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ABSTRACT 

This  is  a  theoretical  study  to  determine  a  practical  design 
formula  for  the  buckling  load  of  simply  supported  rectangular  plates 
under  combined  axial  load  and  normal  pressure. 

The  theory  and  method  used  was  based  on  Samuel  Levy's 
integrated  solution  of  the  non-linear  differential  equations  as  de- 
rived by  von  Karman.     The  length-to-width  ratio  (a/b)  of  the  plates 
considered  was  extended  up  to  a  value  of  4:1.     For  each  length-to- 
width  ratio  a  set  range  of  normal  pressures  was  used  in  obtaining 
the  different  buckling  loads. 

The  study  confirmed  the  conclusion  arrived  at  by  Levy;  that 
is,  the  normal  pressure  increases  the  buckling  load.     It  was  further 
observed  and  concluded  that  for  low  values  of  a/b  (less  than  1 .  0)  the 
buckling  loads  obtained  in  the  manner  described  in  Chapter  III  are 
impractical  to  attain.     Yield  failure  will  probably  occur  before  the 
calculated  critical  load  is  reached. 

In  general,  however,  the  authors'  recommended  buckling  load 
formula  for  a  steel  plate  is  of  the  form, 

cr     ■     nr  .  tv 

c       12(1-//)        b 

where  the  value  of  K  is  scaled  from  Fig.  XXXV.     The  basic  limitation 
on  the  accuracy  of  this  value  of  K  is  the  approximated  expressions )  of 
the  deflection  equation.    The  computed  value  of  ac  must  be  compared 
with  the  yield  strength  of  the  mate  ial  since  the  validity  of  the  theory 
only  holds  true  within  the  elastic  region. 

Thesis  Supervisor:      J.   Harvey  Evans 

Title:  Professor  of  Naval  Architecture 
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I.      INTRODUCTION 


Background 

The  discussion  of  buckling  of  plates  is  usually  based  upon  a 
linear  differential  equation  derived  under  the  assumption  that  the 
deflections  of  the  plates  are  small  in  comparison  with  their  thickness. 
In  particular,  present  design  practice  on  hull  plating  under  combined 
normal  pressure  and  axial  load  is  often  based  on  approximate  critical 
buckling  load  formulas.     The  normal  load  has  generally  been  relegated 
to  minor  roles.     The  reasons  purported  to  supported  this  were  very  well 
founded,  both  from  the  standpoints  of  theory  and  practice.     Samuel 
Levy  and  others  [  2) ,   [3]  arrived  at  the  conclusion  that  normal 
pressure  increases  the  buckling  load  and  thereby  indicated  that  it 
would  be  conservative  design  to  neglect  the  effect  of  normal  pressure. 
The  same  conclusion  was  stated  by  Bleich  [  1] .     Bleich  further-  made 
mention  of  the  fact  that  deflection  of  ship's  hull  plating  usually  does  not 
exceed  one-half  the  plate's  thickness  and  therefore  concluded  that 
linearized  theory  could  be  applied. 

It  is  worthy  to  note  at  this  point,  however,  that  deflections 
exceeding  half  the  plate's  thickness  are  possible,   especially  in  the 
case  of  simply  supported  plates.     In  those  cases  linear  theory  of 
plates  no  longer  applies.    As  will  be  briefly  explained  in  the  sub- 
sequent paragraph,  the  problem  becomes  a  non-linear  stress  problem 
since  deflections  of  the  order  of  magnitude  of  the  plate  thickness  must 
be  considered. 

Briefly  stated,  the  nonlinearity  of  the  differential  equations  has 
its  origin  in  the  fact  that,  in  the  case  of  large  deflections,  there  is  an 
interaction  between  the  membrane  stresses  and  the  curvature  of  the 
plate.    Thi3  interaction  leads  to  non-linear  te^rns  in  the  equations  of 
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equilibrium  of  the  plate  elements. 

The  complete  differential  equations  of  the  problem  were 
formulated  by  von  Karraan,  who  added  those  non-linear  terms 
pertaining  to  the  flexural  rigidity  of  the  plate.     Timoshenko  [  6]  t 
Marguerre  and  Trefftz  [  8 J  derived  the  expressions  for  the  strain 
energy  of  plates  with  large  deflections. 

An  attempt  to  develop  the  large  deflection  theory  for  rec- 
tangular plates  under  combined  bending  and  longitudinal  compression 
was  made  by  Bengston  [  7] .     His  studios  included  plates  with  simply 
supported  and  clamped  edges.     However,   Bengston  introduced,   in 
the  course  of  his  analysis,   certain  arbitrary  assumptions  which  in 
part  contradict  each  othe  r  and  it  is  very  doubtful  whether  the  results 
of  his  analysis  can  be  considered  as  entirely  correct. 

In  a  series  of  papers  Samuel  Levy  gave  exact  theoretical 
solutions  for  rectangular  plates  with  large  deflections.    The 
theories  developed  include  simply  supported  and  clamped  plates. 

The  problem  of  rectangular  plates  carrying  longitudinal 
compression  and  normal  pressure  is  of  prime  importance  in  the 
design  of  the  hull  plating  of  ships.     It  is  closely  related  to  the 
question  of  buckling  strength  of  plates. 

Statement  of  the  Problem 

It  is  common  knowledge  that  factors  of  safety  are  actually 
factors  of  ignorance  which  in  some  cases  have  remained  unchanged 
for  a  period  of  time  for  no  obvious  reason.     Factors  of  safety  may 
be  used  to  account  for  one  or  more  of  the  following:    (a)  material 
imperfections,   (b)  faulty  construction  practices,  (c)  material  de- 
gradations due  to  corrosion,   and  ( i)  lack  of  ability  to  make  rigorous 
calculations.     The  conservativeness  of  design  as  previously  mentioned 
prompted  the  authors  to  ask  in  what  way  would  conservative  design 
affect  the  factor  of  safety;  that  is,   does  conservative  design  actually 
compound  the  factor  of  safety?    If,  for  instance  the  safety  factor  is 
increased  by  only  a  few  tenths  of  a  percent,  then  one  is  justified  in 
using  linear  theory.     However,  if  the  increase  becomes  consiie  able, 
the  design  may  become  uneconomical. 
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With  the  salient  points  stated  in  the  previous  paragraph  in 
mini  and  also  with  the  conclusion  of  Levy  on  the  effect  of  normal 
pressure  on  the  buckling  load  of  plates,  the  authors  came  to  the 
conclusion  that  there  is  a  need  for  further  investigations.     First, 
these  investigations  afford  the  possibility  of  checking  the  accuracy 
of  simpler  approximate  methods  and  formulas.     Secondly,  and 
actually  the  ultimate  goal  of  this  thesis,  these  investigations  would 
hopefully  lead  to  a  formulation  of  a  simple  design  formula  for  plates 
subjected  to  the  simultaneous  action  of  normal  and  edge  loads. 

A  literature  survey  disclosed  there  is  no  simple  formula  which 
a  Jesigner  can  use  with  great  facility.     Fortunately,   general  solutions 
to  the  differential  equations  are  available  in  the  form  of  Fourier  series. 
Though  the  analysis  of  these  solutions  is  highly  involved  and  the  numeri- 
cal work  for  obtaining  special  solutions  is  very  laborious,  a  method 
(Chapter  III)  has  been  devised  which  lends  itself  to  computer  programming. 

This  thesis  is  concerned  mainly  with  plates  simply  supported 
on  all  edges  and  subjected  to  combined  loadings  of  uniform  normal 
pressure  and  axial  load  in  one  direction.     It  is  further  hoped  that  it 
will  serve  as  a  forerunner  of  later  investigations  based  oh  different 
boundary  conditions. 
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II.      THEORY 


Nomenclature    (See  Fig.   I): 

a  plate  length  in  the  x-direction. 

b  plate  width  in  the  y- direction'. 

a    =  a/b 

h  plate  thickness. 

w  plate  deflection. 

x,  y  coordinate  axes  with  origin  at  corner  of  plate. 

E  Young's  modulus. 

iu  Poisson's  ratio.     (Note:    Tabulated  results  and  figures  are 

for  steel  with  11  =  0.  300) 

p  uniform  normal  pressure  on  plate. 

e  average  compressive  strain  at  edges  y  =  0,   b. 

P  axial  load  on  plate. 

F  stress  function. 

D  =  Eh3/12  (1  -  iu2),   Flexural  rigidity  of  the  plate. 
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Fig.   I     Schematic  diagram  of  plate  under  axial  load 
and  normal  pressure. 


Fundamental  Equations 

The  fundamental  partial  differential  equations  governing  the 
deformations  of  flat  plates  have  been  derived  by  von  Karman.    They 
are  given  by  Samuel  Levy  J  2]  and  Timoshenko  ;  6j .    The  exact  math- 
ematical analysis  to  detemine  the  buckling  strength  of  a  simply  sup- 
ported flat  plate  under  combined  edge  compression  and  normal  loading 
involves  the  integration  of  these  equations;  viz. , 

2 


,4„  „4„  .4, 


Tw's 


+  2  — * — 5  +  —-j  *  E 
Ox  9x*  9y*       3y 


.2         '        ,2..      „2. 


J  Ox  y 


(1) 


94w  ,    9     04w      .  84w      p       h  r  02F  02w  .    P>2F  d2w     0   02F      >2w^      /Q. 

rz  +  2  „  2,  2 +  ri ■  b + t>  TJ7J    rr  rr  * 2  0^07  ix-^     (2) 

Ox  Ox  wj        Oy  )y     Jx         Ox     Oy  j  j 

The  boundary  conditions  for  a  simply  supported  rectangular 
plate  to  be  satisfied  by  equations  (1)  and  (2)  are  for  deflection,  w, 
anl  edge  bending  moments  per  unit  length  to  be  zero  at  the  edges 
of  the  plate;  viz. , 

»Jk  ,2     > 

m    ■  -  D  :    ^-?   ♦  t*  --?       ~   0       at    x  *  0,  x  ■  a  (3) 

Ox  yy 

2  2 

m^  tt  _  D      — w   +  u  — k    J    *    0       at    y  «  0,  y  »  b  (i) 

y  Oy*  Ox* 

The  complete  solution  is  more  fully  explained  in  reference  [  2] 
and,  without  loss  of  gene  ality,  the  deflection  equation  is  here  only 
reproduced  in  the  following  Fourier  series, 

-   t     t9m.m-*P—aP  <5> 

m=rl   n»  1 

where  the  undetermined  constants,  w  must  satisfy  the  relation(s) 

expressed  by  equation  (9)  of  the  same  reference. 

The  average  compressive  stain,   e,  at  the  edges  y  »  0,  b  is 
computed  from  equation  (U)  of  [  2]  as: 

or  or 
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In  this  thesis,  the  deflection  equation  for  large  values  of  a/b 
(a  >     1)  is  approximated  by  the  following  expression, 

m  m  w,    .  sin  --  sin  Vf  4  w*  -  sin  3  —  sin  V 
1,1  a  b  3,  1  a  b 

♦  wr    .  sin  ^  sin  5*  ♦  Wf    ,  sin  ?I*  sin  ^    ,       (7) 

o,  1  a  b  <-,  1  a  b 

and  for  a/b  -  values  in  the  vicinity  of  1  and  less  than  1,  the  approximated 
deflection  equation  is, 

••»i,i 8in "»- 8io  £  +  wi.  3 8ln  " sln  tF 

*  w,   ,  sin  !?2  8ln  ^  -  w         8ln  ?JL*  8ln  ?gi 

3,  1  a  b  3, 3  a  b 

#  w,    ,  sin  *£  sin  5fZ  +  w_    ,  flin  ^  sin  VL  .       (8) 

1,5  l  b  5,  1  a  b 

There  are  no  special  reasons  for  using  two  deflection  equations 
for  different  ranges  of  a  values,   except  that  eouation  (7)  restricts  the 
deflected  shape  of  the  plate  to  one  sine  wave  across  the  width  and  a 
combination  of  four  sine  waves  along  its  length.     The  errors  involved 
by  this  approximation  are  expected  to  be  less  than  five  percent  [  2] , 
(  3} .     For  this  reason,   equation  (7)  could  be  reduced  to  a  fewer 
number  of  terms  as  a  is  reduced.     This  is  evidenced  by  the  deflection 
equation  used  fo    or  ■  3.0  in  reference  [  2].     More  specifically,  the 
contribution  of  the  higher-ordered  deflection  coefficients  in  the  x- 
direction  becomes  less  significant  as  a  is  reduced  from  an  at  >~>  1. 

However,  as  a  is  further  reduced  the  assumption  of  one  sine 
wave  across  the  width  of  the  plate  becomes  incorrect  since  the 
higher-ordered  deflection  coefficients  in  the  y- direction  become 
increasingly  significant  and  cannot  be  neglected  unless  the  desired 
degree  of  accuracy  is  sacrificed.     Thus,  for  such  a  values,   equation 
(8)  is  used.     As  before,  it  is  believed  that  the  errors  incurred  by 
using  a  finite  number  of  deflection  coefficients  v/ould  give  results 
within  the  desired  accuracy  (see  Table  1  3  of  [  2] ). 
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The  determination  of  the  value  of  a  applicable  to  equation  (7) 
or  (8)  and  the  more  general  case  of  expanding  the  work  of  Levy  et  al 
on  simply  supported  rectangular  plates  within  the  practical  values  of 
a/b  are  problems  examined  in  this  thesis. 


HI.      PROCEDURE 


With  some  additional  steps  and  explanatory  comments  incorporated, 
the  procedure  followed  in  solving  for  the  deflection  equations  was  essentially 
the  same  as  that  outlined  in  reference  [  3] .     The  steps  were  as  follows: 

1 .  The  family  of  four  or  six  simultaneous  cubic  equations 
corresponding  to  the  same  number  of  unknown  deflection  coefficients 
was  first  obtained.    This  was  done  by  solving  for  the  coefficients  b       , 

P*  H 

defined  by  equation  (8)  of  reference  [  2) ,  and  substituting  them  into 
equation  (9)  of  the  same  reference.     The  results  are  shown  in  Appendix  A. 

2.  Each  of  the  resulting  equations  was  divided  by  h  .  This  step 
nondimensionalized  the  family  of  equations. 

3.  V  lues  of  w,    ,/h,   w,    „/h,  wQ   , /h,   etc.  were  estimated 
corresponding  to  chosen  values  of  Pb/Eh3  and  pb4/Eh4  ••  the  non- 
dimensional  axial  load  and  normal  pressure,   respectively.    This  was 
the  most  delicate  step  since  at  higher  values  of  Pb/Eh3  (normal  pressure 
held  constant)  more  than  one  solution  was  found  possible  (see  [  3]  also). 
In  other  words,  the  estimated  w's,  when  not  properly  chosen,  could  lead 
to  solutions  other  than  those  corresponding  to  a  continuous  change  in 
buckled  form  from  zero  axial  load  to  the  buckling  load. 

It  was  therefore  decided  that  a  step-by-step  method  of  estimating 
the  w's  be  employed.  In  conjunction  with  the  subsequent  steps  4,  5,  and 
6,  this  method  proceeds  as  follows: 

a).      Corresponding  to  three  low  values  of  Pb/Eh  ,  including 

Pb/Eh    =*  0,  three  sets  of  solutions  were  obtained  using  a  reasonable 

3 
set  of  w-estimates.     These  values  of  Pb/Eh    were  hoped  to  lead 

to  unique  solutions.     For  these  estimated  w's,  refe  ence    2]  or 

|  3}  was  used  as  a  guide. 
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b).      From  (a)  above,  succeeding  estimates  of  ws  as 

axial  load  is  increased,  were  based  on  the  continued  trend 

3 
of  each  w-vs-Pb/Eh    curve;  viz. ,  the  slope  and  the  rate  of 

change  of  slope  computed  from  the  last  three  correct  points 

were  employed  in  an  abridged  Taylor  series  of  the  form, 

h  h2 

W(XQ  +  h)  »    W(XQ)  ♦  W  '  (xQ)   jj  +    w*  (xQ)  yi-   , 

where,  w(xQ  ♦  h)  «  next  estimate  of  the  particular  w. 
w(xQ)  a  last  computed  value  of  w. 

w  '  (xQ)       a  slope  obtained  from  the  last  two  computed 


points. 

w*  (xQ)       a  rate  of  change  of  slope  based  on  the  last 
three  computed  values  of  w. 

4.  The  resulting  cubic  equations  shown  in  Appendix  A  were 
linearized  by  expanding  the  right-hand  side  of  each  equation  in  a  Taylor 
series  in  the  neighborhood  of  the  estimated  values  of  w.      /h,  w1    „/h, 
w«   . /h,  etc,  omitting  higher  ordered  terms. 

5.  Crout's  method  [  5]  was  used  in  solving  for  the  difference 
between  the  estimated  w's  and  their  improved  values. 

6.  Step  5  was  repeated  until  such  time  tint  the  calculated 
error  sum  was  less  than  a  test  constant.     The  error  sum  was  defined 

by  the  authors  as  the  measure  of  closeness  of  successive  approximations. 

In  equation  form,  it  is  as  follows; 

N 

Error  «        /       x. 

i=l 

where,    x.    represents  the  absolute  value  of  the  difference  between  an 
estimate  w  and  its  improved  value. 

The  test  constant  used  for  equation  (7)  was  0.0001  and  for 
equation  (8),  0.000001.     Why  two  different  values  of  the  test  constant 
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were  used  would  become  evident  from  a  atu  iy  of  the  results  in 
Chapter  IV. 

7.  The  ave  age  strain  corresponding  to  each  axial  lc  a  i  was 
computed  from  equation  (6). 

8.  In  the  vicinity  of  the  buckling  load  where  the  strain  startel 
to  change  more  rapidly  than  the  axial  load,  the  w1    .  !h  was  made  the 
independent  variable  in  place  of  Pb/Eh   .     The  other  deflection  coef- 
ficients were  computed  in  a  similar  fashion  as  before. 

4         4 
Values  of  the  buckling  load  wee  obtained  for  values  of  pb    /Fh 

equal  to  2. 50,   7.  50,    12. 50,    18. 00,   24.  50,   and  30. 00,   and  at  values  of 

a/b  up  to  4.00.     These  we»e  plotted  versus  a/b  and  analyzed  for  a 

possible  comparison  with  Bryan's  classical  solution. 

The  whole  procedure(s)  had  been  programmed  with  the  use  of 
the  IBM-7094.     Details  not  oth   rwise  covered  here  may  be  founJ  in 
the  programs  shown  in  Appendix  B. 
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IV.      PESULTS 


Tables  I  to  XII  are  sample  results  of  the  computer  programs. 
Tables  I  to  VI  were  obtained  using  equation  (8)  for  an  a/b  =  1.00  and, 
the  rest  of  the  tables  were  for  an  a/b  ■  4.00  using  equation  (7). 

Figures  II  to  XXXIII  are  the  graphical  results  of  the  same 
programs.     The  curves  drawn  are  the  stress-strain  curve  analogy 

for  plates  with  the  non-dimensional  axial  load,   Pb/Eh  ,  as  ordinate 

2     2 
and  the  non-dimensional  average  strain,   eb   /h  ,   as  abscissa.     (Note: 

Not  all  results  for  the  chosen  a/b-values  used  were  plotted. ) 

Tabulated  values  of  the  critical  loads,   rounded  to  the  nearest 
hundredth,  are  shown  in  Tables  XII  and  XIV. 

Figure  XXXIV  is  a  plot  of  the  critical  Pb/Eh3  versus  a/b.     The 
curves  between  zero  a/b  and  the  smallest  value  of  a/b  used  have  been 
extrapolated.     The  curves  have  been  faired  at  points  where  the  results 
of  both  equations  (7)  and  (8)  either  become  tangent  or  intersect.     Re- 
finements of  the  curves  in  the  vicinity  of  the  cusps  were  made  by  ob- 
taining more  results  at  the  a/b-values  concerned.     These  results  are 
shown  in  Tables  XIII  and  XIV. 
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Table  VII    -    Values  of  deflection  coefficients  for  various  values  of  axial 

compressive  load  in  the  x- direction,   P,  for  simply  supported 
rectangular  plate,  a/b  *  4.00,  u  «  0.300.     Normal  pressure, 
P 


2.50  Eh4/b4. 
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Table  VIII     -       Values  of  deflection  coefficients  for  various  values  of  axial 

compressive  load  in  the  x-direction,  P,  for  simply  supported 
rectangular  plate,  a/b  ■  4.00,  \s  =  0.300.  Normal  pressure, 
p  =  7.50  Eh4/b4. 


Pb 

Wl.l 

w3,  1 

W5,l 

W7,l 

eb2 

£? 

h 

h 

.   h 

h 

?~ 

0. 

0.961 

0.223 

0.064 

0.018 

0.  115 

0.  500 

0.986 

0.242 

0.073 

0.021 

0.627 

1.000 

1.014 

0.264 

0.085 

0.024 

1.  143 

1.500 

1.044 

0.288 

0.099 

0.029 

1.663 

2.000 

1.076 

0.314 
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2.  188 
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Table  IX    -     Values  of  deflection  coefficients  for  various  values  of 
axial  compressive  load  in  the  x- direction,   P,  for  sir 

supported  rectangular  plate,  a/b  *  4.  00,  \x  m  0.  300. 

4     4 
Normal  pressure,  p  •  12.  50  Eh  /b  . 


Pb 

Wl,l 

%J> 

W&.1 

W7.1 

6b2 

Eh3 

h 

h 

h 

i 

h2 

0. 
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.037 

.263 

.500 
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T  'ble  X    -       Values  of  deflection  coefficients  for  various  values  of  axial 

compressive  load  in  the  x- direction,   P,  for  simply  supported 
rectangular  plate,  a/b  =  4.00,   <x  =0.300.     Normal  pressure, 
p  »  18.00  Eh4/b4. 


Pb 

Wl.l 

w3,  1 

W5,l 

w7,  1 

«b2 

e7 

h 

h 

h 

-if— 

h2 

0.. 

1.723 

0.464 

0.171 

0.059 

0.448 

0.500 

1.755 

0.485 

0.186 

0.066 

0.984 

1.000 

1.789 

0.507 

0.202 

0.074 

1.524 

1.500 

1.824 

0.529 

0.219 

0.082 

2.068 

2.000 

1.861 

0.552 

0.237 

0.092 

2.618 

2.500 

1.899 

0.575 

0.256 

0.103 

3.173 

3.000 

1.938 

0.597 

0.276 

0.115 

3.733 

3.500 

1.978 

0.620 

0.296 

0.129 

4.300 

4.000 

2.018 

0.643 

0.318 

0.144 

4.873 

4.500 

2.059 

0.665 

0.340 

0.160 

5.452 

5.000 

2.101 

0.687 

0.362 

0.178 

6.039 

10.247 

2.305 

0.795 

0.544 

0.578 

12.930 

10.254 

2.295 

0.792 

0.543 

0.585 

12.956 

10.258 

2.285 

0.789 

0.541 

0.591 

12.977 

10.260 

2.275 

0.785 

0.539 

0.597 

12.995 

10.258 

2.265 

0.782 

0.536 

0.603 

13.009 
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Table  XI  -        Values  of  deflection  coefficients  for  various  values  of  axial 

compressive  load  in  the  x- direction,    P,   for  simply  supported 

Normal  pressure. 


i  ectangular  plate,   a/b«4.00,    u  =  0.300 

p    a 


24.  50  Eh4/b4. 


Pb 
Eh3 


0. 

0.500 

1.000 

1.  500 
2.000 

2.  500 
3.000 
3.500 
4.000 
4.500 
5.000 


w 


1,1 


2.053 
2.085 
2.118 
2.  153 
2.188 
2.224 
2.261 
2.298 
2.336 
2.374 
2.412 


w 


3,1 


\v 


5   1 


w 


0.569 
0.589 
0.610 
0.630 
0.651 
0.671 
0.692 
0.713 
0.733 
0.753 
0.773 


h 

0. 

224 

0. 

238 

J. 

254 

0. 

270 

0. 

287 

0, 

304 

D, 

322 

0. 

341 

c. 

360 

0, 

379 

D. 

398 

7,1 


0.083 
0.091 
0.099 
0.109 
0.119 
0.130 
0.  142 
0.155 
0.169 
0.  184 
0.199 


0.672 
1.216 
1.765 
2.318 
2.875 

4.0C3 
4.  5/4 
6.  150 
5.732 
6.319 


11.964 

2.675 

0.916 

0.607 

0.659 

15.448 

11.973 

2.665 

0.913 

0.605 

0.665 

15.478 

11.978 

2.655 

0.910 

0.603 

0.672 

15. 503 

11.981 

2.645 

0.907 

0.601 

0.678 

15.  523 

11.981 

2.635 

0.904 

0.599 

0.683 

15. 540 

q  lb 
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Table  XII  -      Values  of  deflection  coefficients  for  various  values  of  axial 

compressive  load  in  the  x- direction,  P,  for  simply  supported 
rectangular  plate,  a/b  ■  4.00,  u  ■  0.300.  Normal  pressure, 
p  *  30.00  Eh4/b4. 


Pb 

Wl,l 

w3,  1 

W5,l 

wv 

eb2 

Eh3 

h 

h 

h 

h 

h2 

0. 

2.286 

0.643 

0.261 

0.  101 

0.860 

0.500 

2.318 

0.662 

0.275 

0.  110 

1.409 

1.000 

2.350 

0.681 

0.290 

0.  118 

1.963 

t.  500 

2.384 

0.701 

0.  306 

0.  128 

2.  520 

2.000 

2.418 

0.720 

0.322 

0.138 

3.082 

2.  500 

2.452 

0.739 

0.338 

0.  149  . 

3.647 

3.000 

2.487 

0.759 

0.355 

0.161 

4.217 

3.500 

2.522 

0.778 

0.372 

0.  173 

4.791 

4.000 

2.558 

0.797 

0.39C 

0.  186 

5.  369 

4.500 

2.594 

0.816 

0.407 

0.200 

5.952 

5.000 

2.630 

0.835 

0.425 

0.215 

6.  539 

13.331 

2.920 

0.996 

0.651 

0.724 

17.476 

13.337 

2.910 

0.993 

0.649 

0.730 

17.502 

13.339 

2.900 

0.990 

0.648 

0.736 

17.524 

13.339 

2.890 

0.987 

0.646 

0.742 

17. 542 

13.337 

2.880 

0.984 

0.643 

0.74? 

17. 556 
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Table  XIII    - 


Values  of  (Pb/Eh  )  critical  at  various  values  of  a/b 
and  pb   /Eh  .    These  are  results  using  equation  (8). 


pb4/Eh4  - 

a/b 

2.50 

7.50 

12.50 

18.00 

24.50 

30.00 

0.60 

79.51 

79.60 

79.68 

79.78 

79.89 

79.98 

1.00 

27.56 

28.07 

28.59 

29.17 

29.86 

30.46 

1.20 

16.89 

17.49 

18.06 

18.66 

19.32 

19.86 

1.40 

11.33 

12.29 

13.13 

13.95 

14.84 

15.54 

1.60 

7.90 

9.35 

10.48 

11.57 

12.74 

13.67 

1.70 

6.77 

8.41 

9.68 

10.93 

12.30 

13.43 

1.80 

5.93 

7.71 

9.14 

10.60 

12.31 

13.86 

1.90 

5.36 

7.19 

8.82 

10.62 

13.20 

19.35 

2.00 

4.96 

6.83 

8.72 

11.38 

17.03 

17.65 

2.10 

4.89 

6.54 

8.96 

14.79 

15.68 

16.38 
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Table  XIV    -       Values  of  (Pb/Eh3)  at  various  values  of  a/b 

a         a  critical 

and  pb    /Eh   .     These  are  results  using  equation  (7). 
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4     4 

pb^/Eh 

a/b 

2.  50 

7.  50 

12.  50 

18.00 

24.  50 

30.00 

4.00 

4.20 

6.32 

8.  68 

10.26 

11.98 

13.  34 

3.95 

4.  16 

6.78 

8.77 

10.33 

12.03 

13.37 

3.90 

4.  16 

6.72 

8.87 

10.41 

12.08 

13.40 

3.85 

4.20 

6.67 

8.97 

10.50 

2.  14 

13.45 

3.80 

4.30 

6.61 

9.08 

10.58 

12.21 

13.50 

3.60 

4.  12 

6.49 

9.  64 

11.04 

12.  57 

13.77 

3.  50 

4.03 

6.49 

10.00 

11.32 

12.80 

13.96 

3.40 

3.98 

6.  53 

10.00 

11.65 

13.06 

14.  19 

3.20 

3.92 

6.71 

9.  13 

12.  50 

13.74 

14.77 

3.  10 

3.92 

6.  86 

8.98 

13.09 

14.  17 

15.  13 

3.08 

3.92 

6.89 

8.97 

12.  17 

14.26 

15.21 

3.04 

3.92 

6.97 

8.96 

11.  69 

14.47 

15.38 

3.00 

3.91 

7.05 

8.96 

11.42 

14.70 

15.  57 

2.95 

3.91 

7.  17 

8.98 

11.  19 

15.03 

15.82 

2.90 

3.92 

7.30 

9.03 

11.05 

13.98 

16.  11 

2.85 

3.94 

7.44 

9.  10 

10.97 

13.40 

16.47 

2.80 

3.99 

7.  60 

9.  19 

10.93 

13.08 

15.  19 

2.60 

4.02 

8.44 

9.  79 

11.  17 

12.71 

13.98 

2.  57 

4.04 

6.  89 

9.91 

11.25 

12.74 

13.95 

2.50 

4.  10 

6.  50 

10.  24 

11.  50 

12.87 

13.97 

2.40 

4.20 

6.  33 

10.  83 

11.97 

13.20 

14.  18 

2.30 

4.30 

6.31 

11.  58 

12.60 

13.72 

14.  61 

2.25 

4.36 

6.  33 

9.  15 

12.99 

14.06 

14.90 

2.20 

4.48 

6.  39 

8.81 

13.44 

14.45 

15.26 

2.  15 

4.53 

6.47 

8.64 

11.  64 

14.91 

15.67 

2.10 

4.69 

6.  57 

8.  56 

10.93 

15.43 

16.  16 

2.05 

4.82 

6.70 

8.  56 

10.61 

13.12 

2.00 

4.97 

6.  86 

8.  61 

10.45 

12.  55 

14.30 

1.90 

-- 

7.27 

8.69 

10.45 

12.15 

13.48 

1.70 

6.87 

8.68 

10.06 

11.  36 

12.71 

13.74 
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V.      DISCUSSION  OF  FFSULTS 


The  most  rema   kable  pa^t  of  the  results  w.ia  the  way  the 
critical  buckling  load  behaved  with  respect  to  varying  values  of 
a/b  (Fig.  XXXIV).     A     can  be  seen  in  the  corrolary  figure  (page  64) 
to  Fig.  XXXIV,  there  is  a  striking  resemblance  to  B  yan's 
buckling  loa  1  solution  of  rectanguia     plates  solely  undei  the 
action  of  edge  loadings.     It  may  be  recalled  that  the  wave-like 
pattern  of  B  yan's  curve  of  K  vei  sus  a/b  was  due  to  the  nature  of 
the  K-formula  and  to  the  minimizing  procedure  adopted  in  obtaining 
the  critical  load  for  any  given  value  of  a/b.     Bryan's  K-formula  is 
a  function  of  the  plate's  dimensions  and  of  the  number  of  half-sine 
waves  the  plate  takes  after  buckling.     Thus,  it  would  seem  from  the 
foregoing  statements  to  presuppose  that  the  patterns  obtained  in  this 
thesis  are  also  due  in  part  to  the  number  of  half-sine  waves. 

At  first  glance  the  results  would  give  no  indication  as  to  the 
number  of  waves  the  plate  takes  at  the  incipience  of  buckling.     But, 
checks  had  been  made  by  the  authors  on  the  relationships  between 
deflections  at  the  midwidth  of  the  plate  and  at  various  distances  along 
its  length.     In  other  words  plots  were  made  of  the  deflection  ratio, 
w/h,  versus  the  distance  ratio,  x/a,  at  y  -  b/2.    These  were  done 

for  different  (a/b)'s  at  the  same  normal  pressure.     For  example,  it 

4     4 
was  found  that  for  a  p  =  18.00  Fh   /b     there  were  five  buckles  at 

a/b  =  2.40.     (The  two  half-sine  waves  resulting  from  the  initial 

general  downward  deflection  of  the  plate  due  to  normal  pressure 

were  not  inclu  led  in  the  count. )    When  these  checks  were  expanded  to 

a  wi  ie  range  of  values  of  a/b,  it  was  found  that  points  lying  on  the 

same  trough  of  a  p- curve  gave  the  same  number  of  buckles. 

Since  only  representative  tables  are  given  (Tables  I  to  XII)  it 
was  deemrd  unnecessary  to  insert  all  these  checks  in  this  thesis.     If 
verifications  are  to  be  made  on  this  particular  point  the  authors  suggest 
going  through  the  programs  shown  in  Appendix  B. 
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There  is,  how  eve ' ,  no  positive  explanation  as  to  why  only 
an  o id  number  of  buckles  is  possible.     There  are  either  one,  three, 
five,  etc.,  buckles.     An  even  number  of  buckles,  which  are  obtainable 
in  the  case  of  Bryan's,  seems  entirely  out  of  the  question.     On  this 
particular  aspect,  a  possible  explanation  that  the  authors  believed 
could  have  happened  was  in  the  manner  the  deflection  equations  were 
obtained.     The  approximated  deflection  equations  were  follow-ups  of 
those  use i  in  the  references  \  2] ,   [  3] ,   [  4] .     The  equations  only  have 
the  o  li-numbered,  subscripted  Reflection  coefficients  which,  in  the 
final  analysis,   resulted  in  an  odd  function  for  the  w- equation.    At  this 
point,  there  is  no  way  of  predicting  what  the  results  will  be  if  the  even- 
numbei  ed  subscripts  of  the  deflection  coefficients  are  also  taken  into 
consideration. 

In  reference  to  the  reliability  of  the  results,  it  is  to  be  noted 
that  the  basic  theory  is  vali  i  only  within  the  range  of  Hooke's  law. 
No  allowance  was  made  in  the  original  differential  equations  on  the 
change  of  Young's  modulus  when  the  axial  load  exceeds  the  propor- 
tional limit.     Hence,  all  the  results  of  this  thesis  presupposed  that 
the  material  had  been  within  the  elastic  region  at  all  times.     Experi- 
mental results  are  not  available  on  the  case  of  plates  with  the  3arne 
loading  and  boundary  conditions  such  that,   definite  identifications 
cannot  be  made  as  to  whether  the  results  are  within  the  elastic 
region  or  not. 

When  a  further  gtudy  is  made  of  Fig.  XXXIV  for  a/b  less  than 
or  equal  to  one,  there  is  a  very  wide  discrepancy  between  the  calcu- 
lated critical  loads  and  Bryan's  that  it  is  quite  doubtful  whether  the 
expected  critical  loa is  are  below  the  yield  strength  of  steel  for  a 
majority  of  plate  characteristics. 

The  following  computed  results  are  hoped  to  illustrate  the  point: 

Consider  a  plate  with  a  b/h  =   104.  7  and  acted  on  by  a  pressure, 
p  m   12.    0  Fh4/b4  *  3. 11  psi.     The  b/h  may  correspond  to  a  14.  02- 
pound  plate  with  b  =  36  inches,   and  the  p  may  correspond  to  a  hydro- 
static head  of  7  feet. 
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Case  1  Case  2  Case  3 

a/b-1.50  a/b*1.00  a/b=*0.60 

K»  12.8  K«*  28  K  ■  78 

or    m  31800  psx  a    ■  89000  psi  a    •  192000  pai 

v-  C  C 

If  mild  steel  or  high-tensile  steel  (HTS)  were  used  with  yield 
stresses  of  33000  and  47000  psi,   respectively,  it  is  apparent  that 
cases  2  and  3  exceed  their  corresponding  yield  strengths.     Hence, 
the  theory  fails  for  these  cases.     If  HY80  or  HY100  were  used  then 
theory  is  found  to  fail  only  for  case  3.    At  any  rate,  the  fact  remains 
that  the  basic  theory  does  not  differentiate  between  the  elastic  and 
inelastic  regions.    A  consequent  conclusion,  therefore,  is  the  need 
for  a  check  on  the  yield  strength  whenever  the  formulated  <y    of  this 
thesi  i  is  used. 

Except  for  the  restrictions  mentioned  in  the  previous  paragraph, 
it  is  obvious  that  the  results  obtained  confirmed  the  conclusion  of  Levy 
on  the  effect  of  normal  pressure  on  the  buckling  load.    But,  it  is  also 
evident  that  to  neglect  the  normal  pressure  in  predicting  buckling  strength 
is  unwise  design  from  an  economic  point  of  view. 

Lastly,  the  computer  programs  are  limited  in  scope,  although 
they  worked  properly  in  providing  all  the  needed  results  of  this  thesis. 
They  were  the  consequent  developments  of  the  approximated  deflections 
equations.    A  more  desirable  program  would  be  to  start  with  the 
equations  shown  in  reference  [  2] . 
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VI.      CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

1.  Normal  pressure  always  increases  the  buckling  strength 
of  rectangular  flat  plates.     However,  the  buckling  loads  obtained  by 
means  of  Levy's  solution  are  not  conclusive  values;  that  is,  the 
buckling  load  may  or  may  not  be  within  the  proportional  limit.     In 
other  words,  theory  does  not  differentiate  between  the  elastic  and 
inelastic  regions. 

2.  Plotted  values  of  critical  loads  against  the  ratio,  a/b, 
for  a  constant  value  of  normal  pressure  exhibits  the  same  trend 
and  characteristics  as  the  plot  of  Bryan's  formula  for  zero 
pressure. 

3.  To  neglect  the  effect  of  normal  pressure  in  the  determination 
of  buckling  load  is  unwise  design  from  an  economic  standpoint. 

4.  In  general,  the  buckling  load  formula  can  be  written  in  the 
form, 

*2E        ^2    K 

or       =      —  — r       r-    I        K 
C        120-/)    ^b 

where,  this  time,  K  is  not  only  a  function  of  plate  geometry  but  also 
a  function  of  the  normal  pressure  and  Poisson's  ratio  as  well. 

P  ecommendations 

1.  A  better  representative  plot  of  K  versus  a/b  is  needed  by 
including  the  even-numbered  subscripted   leflection  coefficients  of 
Levy's  solution. 

2.  In  consonance  with  the  fi;st  recommendation,   •  revise  J 
compute*  prog,  am  is  derail  able,   possibly  to  start  right  at  the  series 
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forms  of  the  equations  shown  in  reference  [  2] .     This  step  would 
eliminate  using  approximate  deflection  equations. 

3.     It  is  recommended  that  a  study  be  made  of  the  most 
feasible  experimental  set-up  to  simulate  the  boundary  conditions 
stated  in  this  thesis  and  hence,  the  ensuing  experimental 
verifications. 
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APPENDIX  A 
DFFLFCTION  FQUATIONS 


The  six  (6)  simultaneous  cubic  equations  involved  in  the  solution 
of  the  deflection  equation, 

w       -     w, ,  sin^-sin^+wt_3sm^-si„^- 


+    w3,  I  sta  £**«■  T^  w3.  3  B'«  ^'n  *P 

*  •>,  5  sin  ¥- sin  *F  *  V  1  sin  ^  si"  ?- 


are  as  follows*: 


256  a,      .4  ,  2 

o  * L  £k_  +      4       n  16   Pb 

*6    Eh    iCT^  ***  7  »Tm 

+    <al  W<i-3a2wf ilw13-3a1w^1w3(1 

+    fe a j  +  4a2  +  16a3  +  ag)  wt    j  wj  |  +  (4at+9a2  #  16a3 

*  a4)  wl.  1  Wl  1  »   32a3  Wl,  1  wl.  3  w3,  1  +  81(a4  +  a8)  »l.  1  w3,  3 
'    18(al  *  a8>  wl,  1  %  3  w3,  3  ~  18(a2  +  a4>  wl.  1  w3,  1  w3,  3 

"    (9al  *  64a3  *  25a8)  Wl,  3  W3.  1  "  (?a2  +  64a3  *  25a4>  wi,  3  w3,  , 

♦  (38a1  ♦  36a2  ♦  225a4  ♦  225a8>  w^  3  w^  |f^;  6(a2 

*  3V  Wl,  1  wl,  3  wl.  5  "  6(al  +  3a4>  wl,  1  w3,  1  w5,  1 

♦  162a8  *u  §  w3<  3  wt   5  ♦  162a4  »f   ,  ^  $W$%4  (25^  ♦  4a, 

*J^e  corresponding  significance  of  each  subscripted  a  is  shown  at  the  end  of 
this  Appendix. 


''  '; 


r  i 


or  %  ifi< 

■ 


■ 


0 


I  ,  I  *  b  »J  .  r  ,1 

r.: 


7C 


♦  81  a8  ♦  16  ag>  W%   %  W*   f  ♦  <3a2  ♦  64a3)  w\  %   Wf  , 

+  (3ax  ♦  64a3)  w^  f  W§  ^  (64a3  ♦  274a8>  w%   3  w^  %   w^  § 

♦  (64a3  +  274a4}  wl,  3  w3,  1  W5,  1  "  144  (a3  *  a9>  wl.  3  w3.  3  wl,  5 

-  9(at  ♦  16a3  ♦  49a4  ♦  9a8>  W%   |Wg  gWf  ,  -  9(a2  ♦  16a3 

♦  9a4  +  49a8)w3,lw3,3wl,5*  144(a3  +  a5)w3,tw3,3w5,l 

♦  9(a2  ♦  81a4)  w^  |W^  g  ♦  9<at  ♦  8Ug)  w^  §  w§.  % 

-  (25a0  +  441a.  ♦  256a,.)  w,  ,w2  ,  -  (25a,  +  441a0 

<J  4  0         1,  6      0#  I  1  o, 

♦  256ag)  w^  xw\  $  ♦  729(a4  *  a8>  w3<  gW^w^  |a 

256al    jgb4    .       _4_      hi2  16     Pb  3 

"     3,6       Fh        3<l-,,2)a9Wl'3\2    Ehw1.3-*2wlfl 

♦  (9aj  ♦  4a2  ♦   16a3  ♦  ag)  w2^  1  wt   ,  ♦   16a3  w2^  %  w^  f 

"  9(al  *  V  Wt  1  W3,  3  -  (9a2  *  64a3  +  25a4)  *l.  1  wl  1 

-  2(9a1  ♦  64a3  ♦  25ag)  Wj      Wj    3w3   .  ♦  (Seaj  +   36a2 

♦  225a4  «  225a8)  m%   ,  w3<  jW^i  (8Uj  ♦  a2>  w}#  3 

♦  (272a3  *  625a4  ♦  625ag)  W%   j  w2   ,  -  243ai  w2^  g  w$  3 

♦  9(36aj  •  a2  *  144ag  ♦  9a13>  Wj   g  w2^  g  -  3(a2  ♦  3a8>  w2^  %  «%   % 

♦  2(3a2  ♦  64a3)  wf   ^^w^*  <64a3  *  274a8>  w^  j  w^  f  Wj   5 

♦  (64a3  +   274a4)  Wl.  1  W3.  1  W5,  1  -  144(a3  *  a9>  Wi,  1  w3,  3  Wl,  5 

-  9(at  +   16a3  +  49a4  +  9ag)  wi    xw3   3W5   t  '  <25a2  +  441a4 

♦  256a5)  W%   j  w2   ,  -  (9a2  ♦   64a3  ♦   1225a3)  w2^  jW^ 

"   256a3W3,lW5,l  '  512a3wl,3w3,lWl,5  "  2<9al  +  64a3 
+    1225a4)  wt    3w3    jW5    j  +  9(2a2  ♦  64ag  +  225a4 


8liwI,8w£l!*(8'M8+8,M,+  r  %aw  I  Jw  (  *Ij8)  + 

I  I  o  ,o        £  Si 

nJ  '    <a«38S  )     - 

15  I  I  S 


*<<  -    -m 


3br  *c,xwi  jw(8J  *  £**  +  r*1 


rjwf  i j 


'£  +    ,J3>€ 


r 


:)  ♦ 


)  - 1,. 


8 
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♦  256a9  +  ««tf>*t#  !*!,»%!*  9<2^+  64a34.  256a5 

♦  225a8  ♦    i9a12)  wg   xWg  3w&  |  ♦  (225ax  +  4a2  ♦  256a3 

♦  a10)  wl,  3WL  5  *  (256a3  *  2401a4  *  4096a5 

♦  2401ai2)wl.3wll+  356a3W3,lWl,5w5.1-9(25ftl 

♦  64a3  ♦  9a10)w3<  gW^  g  -  9<«4a3  ♦  144a4  +  169a13)  Wg^w^  &w^  % 

-  (25a  ,  ♦  576a,  ♦  8281a,  „)  w,    _  wK   , . 

256al    pb4        4         h2  144    Pb  3 

♦  16a3  w^  t  wt   3  ♦  (4at  ♦  9a2  *   16a3  ♦  &4>  w2^  jWg^-  9(a2 

♦  a4>  wl.  1  w3,  3  -  (9al  +  64a3  *  25a8)  Wl,  1  wl,  3  '  2(9a2 
+  64a3  ♦  25a4)  wj    jW^^   gWg   -  ♦  (SSaj  +  36a2  ♦  225a4 

4  225a8)  w^  1  w,    3  *3  ,  ♦  (272a3  *  625a4  ♦  625a8>  w2^  3  w^  , 

♦  (aj  +  81^)  wp  j  -  243a2  w2      * i         ♦  9^  +  36a2  +  I44a& 

♦  9a12)w3#1w2<3  ~  3<ai  ♦  3a4)w2   tw51  +  (64a3 

♦  274a8)  w1#  j  wx  ^  (64a3  *  274a,)  w^  ,  w^  8  w5#  § 

♦  2(3at  ♦  64a3)  wt   jWg  jWg  -  «  9(2a2  ♦  16a3  +  9a4 

+  49a8)  Wl,  1  w3,  3  wl.  5  -  144(a3  +  a5>  wl,  lw3,  3  W5,  1 

-  (25at  ♦  441a8  *  256ag)  w^  t  w2  &      256a3  w2^  gw^ 

(9ax  4*  64a3  «•  1225a4)  w2   3w&  j  -  2(9a2  +  64a3 

♦  I225a8)wl3w3  jWl   5>512a3w1   3w3  tw5  1+  9(2a2 

♦  64a3  *  225a4  ♦  256a9  *  49a13>  fjjW,  ,Wj  , 

♦  0(2.1  ♦  64a3  ♦  256a5  ♦  225a8  »  49a12>  w1# ,  W,  ,  w&<  f 
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.j»M£  + 


2    4 


• 


- 

D 


,    t 


U   4 


pflKi  +    r£  (.£   4 


***!    ♦ 


1 4ew 


'V11 


■ 


q    + 
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4    ^B**'0 
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2 
+  (256a3  ♦  2401a8  ♦  4096aQ  ♦  2401at3>  Wg  j  w,  5 

♦  (4a1+  225a2+  256a3+a6)w31w2   1^656a3w1<3w1   5w51 

-  9(25a?  ♦  64a3  ♦  9a0)  w3  gw2   ]  -  9(G4ag  ♦  441ag 

♦  169a12)w3 ,3w1.5W5,l-(25al+   576a3  +  8281a13)  wlt  5  W5,  1* 
256al     pb4  A      324        h2    w  144    Pb  w 

'  ^"  fe   sTTT)  =•   3'3"7"  ^h  3<3 

-  9<*1  +  a8)  wf,  1  wl,  3  "  9(a2  +  a4>  wl,  1  w3, 1  +  31(a4+  a8)wL  I  W3,  3 
+  (38a1  +  36a2  +  225a4^225a8)  wx   1  Wj   gw3       -  81ax  w3  3 

~  81*2W3,  1  *  9(36a!  +»2  +  144a9*9al3)  Wl.  3  w3,  3 

♦  81^  ♦  »2)w3  s  *  9(at  +  36a2  +  144a5  +  9al2>  w|  t  w3  3 

♦  81a8  w2^  t wu%  ♦  81a4  w2^  ,  w.^  §  -  144(a3  ♦  *0)  w^  ,  fj^w,   , 

-  9(2a1*16a3+ 49a4  +  9a3)w1   xwx   gw5   1  -  9(2a2 

+  16aS  +  9*4*49aC)wl,  i  w3#  jw^-  144(a3  +  a5)wu  %  v,^  %  w§  , 

♦  18(a2  +81a4)  w^  §  w^  ^.j  *  I8(a1  ♦  18a8>  w^  ,  w^  $W%  , 
>  729(a4  +  aQ)w1   x  w1   &w5  t  +  9(2a2  +  64a3  +  225a4 

4-  256ag+ 498^)1^^    3Wg   x  w1    &  +  9(2at  +  64a3  +  256a& 

♦  225a8+  49a12)wu  3 w3<  fW§  j  -  9(25ai  *  84a3 

♦  9a1Q)  Wt  rf  5  -  9(25a2-f  64a3^9a6)  w3<  §  w2^  f 

-  9<64a3  +  441a4+  169«13>  *,   ^ ^ W$ j  -  9(64a3 

♦  441a8  ♦  169«12)  w^  |  w^  %m^  %  ♦  81<16as  ♦  81a4 

♦  a8  +  31a10)  W3,  3  WL  5  +  8lU6a3  +  a4  ♦  81a6 


cjsI**  4  „jb  .    :>~   „w  (ft*€  +  ^AtO  4  cja3K)G    - 

I  a. 


TT  ■ 


■ 

c*  ■   ,    jW,    ,  + t,*8/  ♦  tfl£- 


■ 


SI 

W(c     £ 

1  b     A  A 

g#*a>< 
.aI8  4  pboOI8  4  ir.    «w<gv  lit*  4 

.  (QlmlB  4  8«  4 
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+  81»8)w3,3w5.1+2592a3w3,3w5,l  ' 


258al     Pb4  A  4  h2  16     Pb 

0a   ^«  *5T^*rr****~;*ww^ 


-  3(a2  ♦  3a8)  W* f  wf  %  ♦  81a8  ♦  (3a2  ♦  64a3)  w^  ,  w^  % 

♦  (64a3  ♦  274a8)  W%   ,  Wf   3  Wg   ^  9(a2  ♦  81a4)  w^  ,  **   3 
'  H4(a3  ♦  a9)  w^  ,  w^  3  w3  3  -  9(2a2  ♦  18a3  •  9a4 

+  49a8)  wx    x  w^  t  w^  3  -  256a3w^  3^3l'  (9a2  ♦  64a3 
+    1225ag)  wx   3  w2  x  +  9(2a2  ♦  64a3  ♦  225a4  +  256ag 

♦  49a13)  wh  3  w^  t  w3  3  ♦  (25a,  +  4a2  +  81a8  ♦  16ag)  w2^  f  w^  & 

-  2(25ax +441ag4.  256a9)  wx   jWg  Jw1       +  729(a4 

4  a8)  WU  j  w^  3  w5  ,  ♦  (225aT  ♦  4a2  ♦  256a3  ♦  a1Q)  w2^  ,  wf   , 

♦  <256a3  ♦  2401a8  ♦  4096ag  ♦  2401a13>  w2      wt 

4  656a3  w^  3  w3<  jl^j-  18(25*,  ♦  64a3  ♦  9a1Q)  w^  ^  |W|<  , 

-  9(64a3  +441a4  4-  169al3>  w% f *|#  $W^  %  -  9(64a3  ♦  441a8 

♦  169a12)  w^  t  w^  3w^  j  ♦  81(16a3  ♦  81a4  ♦  ag  ♦  81a1Q)w2  3 w^  , 

+  1296a3  W32,  3  w5,  1  '  <25a2  *  576a3  *  8281al2>  wi,  3  W5, 1 

-  2(25a1  ♦  576a3  ♦  8281a13>  w3   i  wi   5  w5   i  +(625aj  +  a2)  w^  & 

+  (1552a3  ♦  28561a  2  +  28561a13>  wl   &  w  2  ., 
n  256al     pb4  4        h2  400    Pb 

-  3<a,  ♦  3a4)w^  t  w,#  f  ♦  81a4  w^  ,  w3>  3  +  (3a,  •  64a3>  w,   ,  w^  , 


:v 


; 

p~£ *  0 


i  .f  ,  :w  i  jw  (8 


i  i*  e  rw  ^s15*251  4 


ij  '  '(8B  + 

l 

x  + 
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♦  <64a3  +274a4)wulWl    3w3    t  ♦  9^  ♦  BU^w^w*  s 

-  9(2aj  ♦   16&3  ♦  49a4  ♦   9&8)  w^  ,  W%  3  Wf  3  -  144(a3 

+  a5>  Wl,  1  w3,  1  w3,  3  "  (9al  *  64a3  *  1225a4>  *1,  3  w3,  1 

-  256a3w1   gw3   t  +  9(2a    ♦  64a3  +  256a&  +  225a8 

♦  49a12) w1(  3 w3>  x wh$  4  (4&1  ♦  25&2  ♦  81a4  ♦  16af))  m\  ,  W%  , 

-  2(25a2^  441a4+ 256a&)w1    t  Wj    3wg   t  +  729(a4 

♦  a8)  wl.  1  W3,  3  wl,  5  +  (256a3  *  2401a4  4-4096a5  +2401a12)w^  $W%  , 

♦  (4.J  *  225a2  ^256a3  ♦  ag)  **  jW^i  656a3  w^  ,  *f  j  w^  § 

-  9(64a3+  441a4  +  lSga^w^  ^  |W|   ,  -  9(64a3  ♦  441a8 

♦  169ai2)  w3,  1  w3,  3W1,  5  "  18(25a2  *  64a3  *  9a6)  w3,  1  w3,  3  W5,  1 

♦  1296a3w2  ^   %  ♦  81<16a3  +  a4  +  81a6  ♦  Blij)  wj  %  W%  , 

-  (25a1+  576a3  + 8281a13)w3   .  w 2       -  2(25a2  ♦  576a3 

♦  8281ai0)w,    ,wt    ,wR   ,  ♦  (1552a0  +  2856la1Q  +  28561a,  Q)w?    _w_    , 

\l       1,  5     \,  b      5#  1  3  14  13      1,  5     5,  1 

♦  Uj  ♦  625a2)  Wj   r 


The  four  (4)  simultaneous  equations  involved  in  the  solution  of 
the  Reflection  equation, 

w  -  wi.  i sin " sin  H- +  w3.  i sin  3F  sin  I1 

♦   wE_  ,  sin  !ga   sin  ^   *  w7<  ,  sin  7-£i   sln  ^  . 


are  as  follows: 


8  gwx       »^ri#)«^  l( 


fw  1WS 


I 


g*I£>  +  *£***--  +  g*£3)8  - 

r(0M-Jll  , t slacks  4  r  tr.<  H   4 

1       EI  I  t  ,1 

.. 

1 


noitoaftel:  sri* 

.   ■  rtfi 


w  -  — 


•- 
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"  3alw'f,lw3,l  +  (ial  +  9a2*16a3  +  a4)wl,lw3,l 

♦  <3al  ♦  64a3)w^  ,  w%  %  -  (3&1  ♦  25^)  w^  %  Wf   f  4  (4&1 

+  25a2+  81a4+ 16a5)w1    t  w2   j  +  (4aj  +  49a2  ♦  256a5 

♦  8U6> wl.  1  w7.  1  *  6<al  *  3a4}  wl,  1  w3, 1  W5,  1  *  2<3al 

♦  64a5)w1   1wf)   jW7   |*  (Saj*   H4a3  ♦  225a4 

♦  9a6)w3,lw5,lW7,r 


0  ■ 


256al    pb4    .  4  h2    m  144    Pb  w  3 

♦  *«*i*»S*1*lWWJ,lW|.l  -  3^  +  3a4>  w2^  ,  w.    t 

3  2 

+  (aj  +  81a2)w3   t  +  (4*^  *  225a2  +  256a3  +  ag)  wg   %  w^  . 

♦  3(ax+  192a3)w2    x  w?    t  +  (4at  ♦  441a2  +  625a4+  16a?)w3   t  w2 

♦  2<3a1  +  64a3)Wl)  x  w3<  (  *ft  %  -  «*,  +  25a4)  *f   %  w3 ,  w^  g 

♦  3(2a.  +  48a„  +  75a.  +  3ac)  w.    ,  wr    ,  w„   ,. 

1  6  4  o        1,10,11,1 


1^   ^¥       3(l-(^    a7 


400    Pb 


w,. 


"J  T!h    w.%  1 


2  2 

-  3(a.  +  3a4)w,       w3   .  +  (4a.  +  25a2  +  81a,  +  lSa^)  wi    i  wk   i 

-  (Saj  ♦  64a5)wt   1  w?    j  +  <3at  +  64a3)  wx      w2       4  v4ax  +  225a2 

♦  256a3+  ag)  w2   xwp    x  +  (at  +  62.c,i*2)w8   -  +  (4a2  •»•  1225a2 

♦  l.296a3  *  a4)  w%  x  w2 ^  +  3(2&1  ♦  48a3  ♦  75a4  ♦  3a6>  w^  x  wf  j  w?<  , 


♦  6(ai^192a3)w31  w5   tw?    r 


r      r  ■      > 


•jw(  4i.twi, 


i  J 


• 


i*  i 


+  ;: 


-  W     ,  ,W(      •  •     +   mfft,    +    . .  . 


III  / 

I 


' 


•         <  .. 


M 


J  t  c 


.. 


* 1,« 


u{         C.P 
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1 1\     *  t< 


I   Sw  '  l^ 
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256al  pb4  .    4     h2  w     784  Pb  w 
^W   3(W^  a~l7  W7,l-"7^  w7.1 

-  (3a-  ♦  64ar)  w2  ,  w,  ,  +  (4a,  ♦  49a0  ♦  256ar  ♦  81afl)  w2   w„  . 

l  1,  1        0,   1  1  Z  b  O  O,  1         7,  1 

-  (3RJ  •  ?.  a^> w1    jWj   t  +  (4ax  +  441a2  t  Ittej  ♦  18a?)  w2    jW?   - 

♦  3(a,  +  192a„)wQ   «  w2       ♦  (4a,  ♦  1225a„  +  1296a0  ♦  a,  J  w2    ,  wn   . 

1  a       5,  1      r-,  I  1  Z  3         14        .,  1       1,  1 

♦  (at  *  2401a2)w^   j  +  3(2at  ♦  48a3  ♦  75a4  +  3a6>  Wj    t  wg   -Wg   .. 


The  subscripted  a's  have  the  following  significance: 

a-  *  a  aQ  =»  a   /(9c    ♦   1) 

2 

a2  ■  1>«2  a1Q  *  a  2/(\Qa2  ♦  1) 

2  2 

a3  *  or2/(tf2  +   1)  an  -  a2!(2la2  +  I) 

2  2 

a4  ■  a2/(*2  ♦  4)  a12  »  a2/(4<*2  +  9) 

2       2  2  II2 

a5  *  o   I  {a    ♦  9)  «13  =a   /(9a    *  4) 

2 
a6  *  o2/(o2  ♦  16)2  aH  -  *2/(*2  t   36) 

2 
a?  *  a2 /(a2  +  25)2  a15  *  a2  fa2  ♦  49) 

a8  »  a2(4a2  ♦   1) 


where  a?  «  a/b. 


It  » 

c 
I  ,c      I  ,8  a 


«    rt.« 
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•  eB 

>>ro 


-a.»  » *  t.* 


( r  - )  o  *    ■ 
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M3888-3401,FMS»DEPUG»1»1»500,0 


BRILLANTES-LIVADAS 


SOLUTION  OF  THE  DEFLECTION  EQUATION 


THIS  METHOD  IS  SET  UP  FOR  A  CONTINUOUS 


°-     b 
ex 


CHANGE  IN  AXIAL  LOAD  cROM  ZERO  TO  A  VALUE  WITHIN 
THE  VICINITY  OF  THE  CRITICAL  BUCKLING  LOAQ.   ZERO 
LATERAL  LOAD  IS  EXCLUDED. 

MAIN  PROGRAM 


XEQ 

LABEL 

L  I  S  T  8 

FORTRAN 

READ  1»   IND 

1  FORMAT  f 151 
PRINT  2   • 

2  FORMAT  (82H1   PNORM    AXIAL 


W  {  1  ,  1  )  /  H   a'  (  3  »  1  )  /  H   W  (  5  ♦  1  )  /  H 


1D/H   STRAIN         ERROR   ) 
DO  44   JOB  =  1 ,:ND 

READ  3»  U*  ALPHA*  PNORM,  LOOP,  LAP 
3  FORMAT  (3F8. 0*213) 

DIMENSION   Bl(4),  B2(4),  B3(4) 

CALL  CONST  (U,  ALPHA*  PNORM,  Al*  A2,  A3,  A4,  A5»  A6,  A7,  A8,  A9, 
1  A10,  All*  M  2 »  A13*  A  1  '+  ,  A15,  A 1 6  »  A 1  7  ,  Bl,  B2»  B  3  ) 

INITIAL  ESTIMATES  OF  .v  (  I  )  i  S 
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DIMENSION  W(4) 

W( 1 )  =  0.400 

W(2)  =  0.300 

W ( 3 )  =  0.200 

W(4)  =  0.100 

THE  rOLLOWlNG  IS  THE  ROUTINE  EOR 
CALCULATING  THE  DEFLECTION  CuEFF  I  C  I  EP.TS  FOR 
CHOSEN  VALUES  OF  AXIAL  LOAD  AND  NORMAL  LOAD 

DIMENSION   A(4,4),  C(4),  X(4)»  Z(3)t  S(3»4) 

AX  =  Q  . 

MG  =  1 

10   LL  =  1  »3 

4  CALL  SESAME 

CALL  CROUT  (A,  C»  X,  M) 
GO  TO  (5,  42) »  M 

5  fdooo  =  o. 

DO  6   J  =  1  ,4 

6  ERROR  =  ERROR  +  ABSF(XU)) 

IF  (ERROR  -  0.0001  )   9,9,7 

7  DO  8   1=1,4 

a  w  ( i )   =  w  ( i )   +  x  1 1 ) 

GO  TO  4 
9  CALL  CPRNT  (PNORM,  AX,  W*  ERROR,  Al) 
CALL  STORE  (AX,  V.,  Z,  S,  LL) 
AX  =  *X  +  0.5 
10  CONTH 

n^L  =  0.5 
KO  =  1 
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->?    GO  TO  (23i  25,  27)  ,  KO 
2?  D^L  *  0.05 

DO  24   I  =  1,4 

W(  I  J  ■  S(2,  I  ) 

24  5(3*  !1  =  W(  I  ) 
AX  =  Z( 2) 
Z(3)  =  AX 

AX  =  AX  +  DEL 

KO  =  2 

KK  =  0 

CALL  FXCH  (7,  S) 

GO  TO  13 

25  DFL  =  n.Ol 

DO  26   I  =  1,4 
W(  I  )  =  S(2,I  ) 
2  6  5 ( 3 ,  I  )  =  W (  I  ) 
AX  =  Z ( 2  ) 
7(3)  =  AX 
AX  =  AX  +  DEL 
KO  =  3 

kk  =  n 

CALL  FXCH  (Z,  S) 
GO  xn  13 


THE  FOLLOWING  IS  THE  ROUTINE  WHEN  IN  ' 
VICINITY  OF  THF  BUCKLING  LOAD.   W(l)  IS  MADE 
THF  INDEPENDENT  VARIABLF  IN  PLACE  OF  THE  AXIAL 
LOAD. 


27  DO  2 


1  =  1,4 
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-  n  , 
DO  32   J  =  1»4 

32  ERROR  =  ERROR  +  APSF(XU)) 

IF  (ERROR  -  0,0001)   36,  36,  33 

33  KK  =  KK  +  1 

IF  (KK  -  LOOP)   34,  34,  40 

34  DO  35   I  =  2,4 

3  5  W (  I )  =  W ( I )  +  X (  I  ) 

AX  =  AX  +  X (  1  ) 

GO  TO  30 
^6  TALL  CPRNT  (PNOPM,  AX,  W,  FRROP,  Al) 

IF  (AX  -  Z (2 ) )   ':-4»  44,  3  7 

37  7  (3 )  =  AX 

DO  38   I  =  1,4 

38  S(3,I)  =  A ( I) 

DL1  =  7(3 )  -  Z (2 ) 

39  CONTINUE 
GO  TO  44 

40  PRINT  41 

41  FORMAT  (40H    INFINITE  LOOP  OR  NEEDS  MORE  ITERATION) 
CALL  CPRNT  (PNORM,  AX,  W»  ERROR,  Al ) 

GO  TO  44 
4  2  PRINT  43 

43  FORMAT  (19H    CROUT  WENT  CRAZY) 

44  CONTIN'  I 
CALL  EXIT 

COMMON  W,  AX,  Bl»  B2»  B3 *  Al,  A2i  A3,  A4,  A5,  A6,  A7,  A8»  A9,  A 10, 
1  AllfA12tA13tAl4»A15»A]  •  \ 1 7 •  I   •  '■  »  C 


END 
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SUBROUTINES 

SUBROUTINE  FOR  SOLVING  THE  LINEARIZED  SYST 
OF  EQUATIONS  RY  Tmr  CROUT  REDOX  METHOD 

LABEL 

LIST8 

FORTRAM 

SUBROUTINE  CROUT  (A,  C»  X,  M) 

DIMENSION   A(4,M,  C(4)»  AA(4,4)»  CC(M>  X(4) 

DO  100   I  =  1,4 
100  A  A (  I  » 1  )  =  A (1,1) 

DO  101   J  =  2,4 

A  A  (  1  ,  J  )  =  A  (  1  ,  J  )  /  A  (  1  ,  1  ) 

DO  K2   I  =  2,4 

DO  102   J  =  2»4 

102  A  A  (  I  »  J  )  =  0  . 
J  =  2 

103  II  =  J 

DO  105   I  *  1 1  »4 

LIM1  =  J  -  1 

DO  104   K  =  1 »  L  I  M  1 

104  AA(I»J)  =  AA(I»J)  +  (AA(I.<)  *  AA(<»J)) 
10?  AA ( I » J)  =  A(I,J)  -  AA( I »J) 

IF  ( A  A ( J  ,  J  )  )   106,  116,  106 
10  6  IF  (4  -  J )        no,  no.  107 
107  T  =  J 


84 


J  =  J  +  1 

JJ  =  J 

DO  109   J  ■  JJ»4 

LIM2  =1-1 
DO  108   K  =  1»LIM2 
1 0  8  A  A  (  T  »  J  )  =  A  A  (  I  t  J  )  +  (  A  A  (  I  »  K  )  *  A  A  (  <  t  J  )  ) 

AA( I » J)  =  A ( I » J)  -  AA( I » J) 

109  AA(T»J)  =  AA( I »J)/AA( I »I ) 
J  =  T  +  1 

GO  TO  103 

110  DO  1 11   1=1,4 
X( I )  =  0. 

111  CC( I )  =    n. 

CCM  )  =  C( 1) /AA(  1  ,1  ) 
DO  113   I  =  2»4 
LIM3  =  T  -  1 
DO  1 1?   K  =  1  »LTM3 
it?  CC(I)  =  CC(I)  +  (AA(I,K)  *  CC(O) 
CC( I  )  =  C(  I  )  -  CC(  I  ) 

113  CC(I J  =  CC(  I  )/AA(  T  ,1  ) 
X(4)  =  CC (4) 

DO  115   I  =  1»3 

IT  =  4  -  I 

LIM4  =  II  +  1 

DO  114   <  =  LIM4»4 

114  XUI)  =  X  (  I  I  )  +  (  A  A  (  I  T  »  K  )  *  X  (  K  )  ) 

115  X{  I  I  )  =  CC<  I  I )  -  X (  I  I  ) 
M  =  1 


GO  TO  118 

116  PRINT  117 

117  FORMAT  (28H    SINGULARITY  -  NO  SOLUTION) 
M  s  2 

118  RETURN 
FND 

SUBROUTINE  FOR  CALCULATING  THF  NON-DIMENSICNAL 
STRAIN  AND  THEN  PRINTING  ALL  THE  NEEDED  RESULTS 

LABEL 

LIST8 

FORTRAN 

SUBROUTINE  CPRNT  (PNORM,  AX,  W,  ER,?OP,  Al) 

DIMFNSION  W(4) 

SUM  =  0. 

DO  200   IJ  =  1,4 

XY  =  2  *  IJ  -  1 

SUM  =  SUM  +  (XY  *  WtIJ))**2 

STRAIN  =  AX  +  3.1416**2  *  SUM/(8.  *  Al) 

PRINT  201,  PNORM,  AX,  (W(I)»  1=1,4),  STRAIN,  FRP^P 
201  FORMAT  (F8.2,  F11.5*  4F10.3,  F11.5,  F12.5) 

RETURN 

SUPROUTINE  FOR  STORAGE  OF  THP"r  SETS  OF  VALUES 
OF  AXIAL  LOAD  AND  THE  FOUR  DEFLECTION  COEFFIC 

LABEL 

LIST8 

F  0  R  T  R  A  N 

SUBROUTINE  STOR^  (AX,  W,  Z,  S,  LL) 
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DIMENSION  W(4)»  Z(3)»  5(3,4) 
Z(LL)  =  AX 
DO  20?   I  =1,4 
203  S(LL,  I  )  =  W(  I  ) 
RETURN 


END 


THIS  SUBROUTINE  IS  JUST  A  SHIFTH      DCEDURE 
OF  VALUES  FROM  ONE  STORAGE  LOCATION  TO  ANOTHER. 


LABEL 
LI  ST  8 
FORTRAN 

SUBROUTINE  EXCH  !Z,  S) 
DIMENSION  Z ( 3 ) ,  S (3 ,4) 
DO  2 [J U         I  =  1,4 
?  0  4   J  =  1  ,  2 
JA  =  J  +  1 

204  S( J,  I  )  =  S(  JA,  I  ) 
DO  205  IJ  =  1,2 
JB  =  IJ  +  1 

205  Z( I  J)  =  Z( JB) 

3  2  06   J  I  *  1  »  4 
i  0  6  5  (  3  ,  J  I  )  =  0 . 
Z(3  )  =  0. 
r  C  TU 
END 


SUBROUTINE  FOR  COMPUTING  ALL  THE  NECESS^ 
CONSTANTS  IN  THC  MAIN  PRO' 


LABEL 
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LI' 

FORT- 
SUBROUTINE    CONST     (U*    ALPHA*    PNORM»    Al»     A2*     A3*    A4 •    A3*  •  ■ 
1A9*    A10*    All*    A12*     A13*    Al4»     A15l     U6«     A 17,     51*    B2»     D 3  ) 
DIMENSION    Bl(4)»     B2C4)i    B3(4)»     D(9) 
P(l)     =     AL^HA**2 
d  (  2  )     =     1  •  /  P  (  1 ) 
DO    n ""        1=1*7 
Yl     a     1**2 
11=1+2 
PHI  )    =    P(  1)/(P(1)    +    Yl)**2 

A  1  =     P  (  1    ) 

A2  =     P( 1  )     +  P(2  ) 

A3  *    PCI)     +  81.     *    PC2  ) 

A4  =P(1)+  62  5.     *    P(2  ) 

A5  =    P( 1  )     +  24(  1.     *     P(2  ) 


A6       ■  3«  »    Pfl 

A  7       =  3  .  *    D  (  1 

A  P         *  3«  *     P { ] 

A9        *  3.  *     Pf  1 

AlO    =  3.  *    P (  1 

AH     =  4  .  *     P  (  1 

A 1  2     =  4  .  *     P  (  1 

A  1  i      =  A  .  #     P  (  ] 

A  1  4     =  4  .  *     P  (  1 

A 1 5    =  4 .  *    P  (  ] 

A 1  5     =  4  .  *    P  (  1 

A 1 7     =  6 .  *     PC 


+  64.     *    P ( 3) 

+  5  76.     *    P(3) 

+  D,      *      D(4) 

+  2  5.     *    ° ( 4 ) 

+  6  4.     *    P ( 5 ) 

+  Q .     *    P  {  2  )  +          16«*H(3)+P(4) 

+  225.    *    P ( 2 )  +       2 5 6 •    *    P( 3  J     +    D ( A ■  ) 

+  1 2 2 5 •    *    P ( 2 )  +    1 2 9 6 •    *    P  (  *  )     +    P ( 8 ) 

+  2  5  .     *    P  (  2  )  +                »*P<4)+16.*P 

+  4  41.     *    D  (  ?  )  +        6 :'  5  .     *     P  (  4  )     +     16.     *    P  (  7  ) 

+  40.     *    P  (?)  +        2  56.     *    P  t  5 )     +     8  1.     *     P ( 6  ) 

+  144,    *    p ( 3 )  +       225.    *    P ( 4 )     +       9 .    *    P(b) 
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)  3  0  ]    1=1,4 
12  «  2  *  I  ♦  3 

Y2  •  2  *  I  *  1 

B1(I)  =  2  56.  *  P(l)  »  PN0RM/fY2  »  3*1416**6) 
B2  U  )  ■  4  .  /  (  3  .  Mli  -  U**2  )  *  P  (  I  2  )  ) 
30]  R3(I)  =  -(16.  *  Y2**2J/( 3*1416**2) 

RCTURN 
END 

SUBROUTINE  FOR  COMPUTING  THE  COEFFICIENTS 
C  FOR  THE  USE  OF  THE  CR  )UT  REDOX  METHOD 

*  LAREL 

*  LI  ST  8 
FORTRAN 

SUBROUTINE  SES* 
COMMON  W»  AX,  Bl»  B2»  B3*  Al,  A2»  A3,  A4,  A5»  A'»  A7,  A  •    ,  AlOi 

]         AM,  A  1 .2  »  A  1  3  »  A  1  4 ,  A  1  3  ,  ATA,  417,  M  G »   '  «  C 

DIMENSION    W(4)»     BK4)»     B214H     B3(4)*     A(4»4)»     C(4) 

GO    TO     ( 401  ♦     40  0)  i     MG 
40n    M1>1  )»B3t  1  J*W(1  ) 

A (2*1  )=R3( 2)*W(2 ) 

A(3»l  ) =R3 ( 3  )*W(3  ) 

A(4»  1  )=B3  (  4)*W  (/,) 

GO    TO    40  2 
401     A(i,i)=R?(l)+BMl)  *AX+3«*A2*W(  ]  )  **2~6«*A1  *  '■•'  (  ]  )*W(  ■>  )+Ai  ]  ]     *2 

i  +A14*W(3)**2  +  A16*W(4)**2-2«*A8  «-W  (  2  I  *W  (  3  )  -2  •*AlO*k  (  3  )  *W  (  u  ) 

a  (  ?  ,  i  )  =--3.  ,*A1  *W(  1  )  **?+?.  *A  ]  i  *W  (]  )*V  (  ~>  )  -■?  ,*A8*W(  ]  )  "  tf(3) 
1  +2 • * A6* W ( 2 ) *  W  ( 3  )  -2  •  "  ft  9  *  W  (  2  )  *  W  (  4 )  +A 1 7*  '■''  (  ^  )  *  W  (  4  ) 

M3»l )«-2 **A8*W ( 1 )*W( 2)+2.*A14*W( 1 1 *W ( 3  ) -/ • *A 10*W ( 1 )*W(4) 
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1         +A6*W(? )**2+A17*W( ?  )*W< 4 ) 
A  (h  >  1  )  =-?#*AlO*W(  1  )  *W(  3  )  +  2  ,*A16*W  (  1  )*W(4  )-Ao#W(  ?  )  -" 
40  2  A( 1 »2 )*-3«*Al*Wf 1)**2+2»*A11*W(1 J*W(2 )  +2.*A6*W(  2)*W(3) 

1         -2.*AQ*W(2)*W(4)-2.*A8*W(1)*W(3)+A17*'a'(3)*W(4) 
A ( 2» 2 ) =B2  <  2 )+B^ ( ? )*AX  +  Ai 1*W( 1 ) **2+3«*A3*W< 2)**2+A12*W( 3 )**2 

1        +A1 5*W( 4  )**2+2»*A6*W< l ) *W( 3  ) -2 »*A9*W ( i  >*W ( 4 ) 
A(3»2)*-A8*W(1 >**2+2«*A6*W( l )*W(2 )+2«*Al2*W(2 )*W(3 ) 

1         +A17*W( 1 »*W ( 4 ) +2.*A7*W ( 3 } *W IL ) 
A  (4*  ^  )  =-?.*AO*W(  1  )*W  (  2  )+2.*Al  R#W(  2  )*W(  4  ) +*  7*W  (  "5  )  *#2+A]  7*W  (  1  )  *W  (  3  ) 
A(  1  »3  )  =  A6*W(  2  )  *#?  +  2  ,*A14*W  (  1)*W(3  )-2.*A8*W  (  1  )  *W  (  ?  ) 

1        -2«*A10*W( 1)*W(4)+A17*W(2 )*W(4) 
A ( 2  »  3 ) =-A8*W ( 1 ) **2  +  2 • * Al 2*W ( 2 ) *W ( 3 ) +2 • *A7^  W ( 3 ) * W ( 4 ) 

1        +  7 ,#A6*W ( 1 ) *W ( ? ) +A1 7*W ( 1 ) *W ( 4 ) 
A(3»3)=B2<3 )+B3 (3)*AX+Al4*W( 1 )**2+Al2*W(2)**2+3«*A4*W(3 )* 

1         +A13#W(4)  Jt*2+2.*A7*W(2  )  *W(4) 
A<  4»^  )  =-Al  0*W(  1  )**2+2«*A7*W(  2  )  *'A'  <  3  >+2«*Al3*W(  3  )  *W  <  4)  +A17*W  (  1  )  *W  (  2  ) 
A  (  i  »4)s-A9*W(2  )**2  +  2«*Al6*W(  ]  )*W(4)--?#*Al0*W(  1  )  *W(  3  )+A17*W  (  2  )  * 
A(2»4)=A7*W(3) #*2+2»*Al5*W( 2)*W(4)-2 »*A9*W ( ] ) *W ( 2 )+Al7*W ( ] ) *W ( 3 ) 
A(3»4)=-AlO*W( 1 )**2+2«*Al3*W( 3)*W(4)+A17*W( 1 )*W(2)+2«*A7*W(2)*W(3 ) 
A(4»4)=B2(4)+B3(4)*AX+A16*W(1 )**2+A15*W(2)**2+A13*W(3 )* 

1        +3,*A5*W(4)**2 
r  ( 1  ) = R 1  ( ]  ) -B2 ( ]  I *W ( ]  ) -P 3 ( i  ) *W ( 1 ) *AX-A 2*W ( 1 ) **3  +  3 , *A1*W ( 1 ) *  ?   ( 2 ) 

i      -A11*W(  1 )  *W  (  7  )  **2~A6*  W  (  ^  )  **2*W  (  ^  )  +A9*W  (  2  )  **2*W  (  4  ) 

?      -A14*W(1 )*W(3 ) **2~A16*W( 1 ) *W (4) **2+2»*A8*W ( ] ) *W ( 2 ) *^  (3 ) 

3      +?.*A10*W (1 )#W( 3 ) *W(4 ) -A17*W (  2  )  *W ( 3 ) *W ( 4 ) 
C( 2)=B1 (2)-B2<2)*W(2 )-B3  <2 )*AX*W( 2 )+A]  *W(  !  )  **3-Al  n  ;- ■'  I  ]  I 

]      +A8*W ( 1 )**2*W (3 )-A3*W( 2)**3~A] 2*W(2) *W ( 3 )**2~A7*W< 3  H      (4) 

7      -A15*W(2)*^'(4)**2-2«*A6*W(l)*W(2)*W(3)+2«  4 ) 
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■x  _A17*W(  1 )*W( 3 ) *W(4) 

C( 3 ) *B1 ( 3 ) -B2  C 3 )*W( 3 ) -B3 ( 3 ) *AX*W ( 3 ) +A8*W ( 1 ) **2*W ( ?  ) 

1  -A14*W ( l ) **2*W ( 3 )+A10*W( 1 ) **2*W< 4 ) -A6*W( l ) *W( ?  ) **? 

?  -A12*W(  ?)  **2*W.(  31  -A^*W(3  )  **3-Al^*W(  3)*W  (  4)  **2 

3  -A17*W ( l ) *W( ?  ) *W ( 4 ) -? »*A7*W ( 2 ) *W (3)*W( 4) 

C  (  4  )  =P1  (4)  -B2  (  4  )*W  (4  )  -B*  U  )*AX*W(  4  )+ft  ]  0*W<  1  )  *#2*W(  3  ) 
i  -A16*W( ] ) **?*W( 4)+A9*W( i >*W( 2)**2-Al5*W( 2)**2*W(4) 

2  -A7*W( 2 )*W(3)**2-A13*W(3l»*2*W<4)-A5*W(4)**3 

3  -A17*W( 1 )*W( 2 )*W( 3 ) 
RETURN 

END 


M3888-3401»FMS»DFBUG»1>1»500>0 


91 


SOLUTION  OF  THE  REFLECTION  EQUATION 

THIS  METHOO  IS  SCT  UP  ^OR  A  CONTINUOUS 
CHANGE  IN  AXIAL  LOAD  FROM  ZrPO  TO  A  VALUE  WITHIN 
THE  VICINITY  OF  THE  CRITICAL  BUCKLING  LOAD.   ZERO 
LATERAL  LOAD  IS  EXCLUDED. 


MAIN  PROGRAM 


XEO 

LABEL 

LIST8 

FORTPAN 

REAP  1»   IND 

1  FORMAT  (15) 
PRINT  2 

2  FORMAT  (114H1   PNORM      AXIAL      W(l»l)/H      U»3)/H     W(3»l)/ 
1H     W(3*3)/H     W(lt5)/H     W(5»l)/H     STRAl  ERROR  ) 

DO  44   JOB  =  1  * IND 

READ  3»  U»  ALPHA*  PNORM*  LOOP*  LAP 

3  FORMAT  (3F8. 0*218) 

DIMENSION  W(6)»  R(54)»  Bl(4)»  B2(5)»  B3(3)*  A (6*6)*  C(6)»  X(6) 
CALL  CONST  (U»  ALPHA*  PNORM*  R»  Bl»  B2»  B3) 
INITIAL  ESTIMATES  OF  W(I)»S 
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'.-'  (  1  )  =  0,100 
l\    =  0.010 

W ( 3 )  =  0.010 

W  (  4  )  =  0  .  0  n  l 

W( 5  )  *  0.001 

W(6  )  =  0.001 

THE  FOLLOWING  IS  THE  ROUTINE  FOR 
CALCULATING  THE  DEFLECTION  COEFFICIENTS  FOR 
CHOSEN  VALUES  OF  AXIAu  LOAD  AND  NORMAL  LOAD 

DIMENSION  Z(3) »  S(3»6) 

AX  e  0. 

MG  =  1 

RONE  =  R(l) 

DO  10   LL  =  1»3 

4  CALL  SESAME 

CALL  TROUT  (A*  C»  X,  M) 
GO  TO  (5,  42  )  ♦  M 

5  ERROR  =  0. 

DO  6   J  =  1  »6 

6  ERROR  =  ERROR  +  ARSF(XU)) 

IF  (ERROR  -  0.000001)   9»  9*  7 

7  DO  8   I  =  1*6 

8  W  (  I  )  =  W  (  I  )  +  X  (  I  ) 
GO  TO  4 

9  CALL  CPRNT  (PNORM*  AX,  W<  FRROP,  RONE) 
CALL  STORE  (AX,  W»  7,  5,  LL) 

AX  =  AX  +  0.5 
10  CONTINUE 
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DFL  =  0.5  ' 
KO  =  1 
KK  =  0 

11  DO  12   I  =  1,6 
Fl  =  (  S  (  2  »  I  )  -  S(lfZ))/(Z(2)  -  Z  (  1  )  ) 
F2  =  (S(3,I)  -  S(2»I))/(Z(3)  -  Z(2)) 
F3  =  (7(3  )  -  ZC))/2, 
F  4  =  (  F  2  -  F 1 ) / F 3 

12  W(I)  =  We  I)  +  F2  *  DEL  +  0.5  *  F4  *  DEL**2 
CALL  EXCH  (7,  S) 

13  CALL  SFSAME 

CALL  CROUT  (A,  C»  X,  M) 
GO  TO  (14,  42) ,  M 

14  ERROR  =  0. 

DO  15   J  =  1  ,6 

15  ERROR  =  ERROR  +  ABSF(XU)) 

IF  'ERROR  -  0.000001)   ]9,  IP,  16 

16  KK  =  KK  +  1 

IF  (KK  -  LOOP)   1 7,  17,  40 

17  DO  IB   1=1,6 

18  W(I)  =  W (  I  )  +  X  (  I  ) 
GO  TO  13 

IP  IF  ( W ( 1 )  -  S ( 2  , 1  )  )   2  2,  22,  20 
.20  CALL  CPRNT  (RNORM»  AX,  W»  ERROR,  RONE) 

2(3)  <=  AX 

no  21   1=1,6 
21  S ( 3  ,  I  )  =  W  (  I  ) 

AX  =  AX  +  DEL 
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22 

23 


24 


25 


26 


KK  =  0 

GO  TO  11 

GO  TO  (23,  25.  27) ,  KO 

DEL  =  0.05 

no  ?U       I  =  1»6 

W(  I  )  =  S (2  ,  I  ) 

S  ( 3  1 1  )  =  W  (  I  ) 

AX  =  1(2) 

Z(3)  =  AX 

AX  =  AX  +  DEL 

KO  =  ? 

KK  =  0 

CALL  EXCH  (Z,  Sj 

GO  TO  13 

DEL  =  O.01 

DO  26   1=1,6 

W  (  I  )  =  S  (  2  ,  I  ) 

S 1 3  »  I  )  =  W  (  I  ) 

AX  =  Z ( 2) 

Z(3)  -    AX 

AX  =  AX  +  DEL 

KO  =  3 

KK  e  0 

CALL  EXCH  (Z,  S) 

GO  TO  1 3 

THE  FOLLOWING  IS  THE  ROUTINE  WHEN  IN  THE 
VICINITY  OF  THE  BUCKLH      *D.   W(l)  IS  MADE 
THE  INDEPENDENT  VARIABLE  IN  PLACE  OF  THE  AXIAL 


»5 


LOAr. 
27  DO  28   I  =  1»6 
S(3,  I  )  =  S(2.I  ) 
S  (  2  ♦  I  )  =  S  (  1  ♦  I  ) 

2  8  W (  I  )  =  S ( 3  ♦  I  ) 

1(3)     =    2(2) 

1(2)     ■  1(1) 

AX  =  7.(3) 

TH^  FOLLOWING  VALUES  0^  DL1  AND  DL2  ARE  ARBITRARY. 
THF  FIRST  VALUE  OF  DL1  IS  SO  CHOSEN  SUCH  THAT  WHEN 
W(l)  IS  DECREASFO  BY  DL2»  THE  CONSEQUENT  ITERATED 
SOLUTION  WILL  NOT  GO  BACK  ON  THE  PREVIOUSLY  COMPUTED 
PATH.    (NOTE  -  THF  W(l)  OR  THE  FIRST  HARMONIC  STARTS 
TO  DECREASE  IN  VALUE  BEFORE  THE  CRITICAL  BUCKLI 
LOAD  IS  REACHED) 

DL1  =  0.20 

DL2  =  0,010 
.  MG  =  2 

THIS  LAST  DO  LOOP  IS  HOPED  TO  GO  BEYOND  THE  CRI- 
TICAL RANGE.  START  WITH  LAP  =  30. 

DO  39   LAST  =  1 »LAP 

KK  =  n 

AX  =  AX  +  DL1 

W( 1 )  =  W( 1 )  -  DL? 

DO  29   J  =  2,6 
2  ^  W  (  I )  =  W  (  I  )  +  ( S  (  3 » J  )  -  S  (  2  *  J  )  ) 

CALL  EXCH  (Z»  S) 
30  CALL  SESAME 
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CALL  CROUT  (A,  C«  X,  Ml 
GO  TO  ( 31  ,  42)  ♦  M 

31  cRROR  =  0. 

D  0  3  2   J  =  1  f  6 

32  ERROR  =  ERROR  +  APSF(XU)) 

IF  (ERROR  -  0.000001)   36i  36,  33 

33  KK  =  KK  +  1 

IF  (KK  -  LOOP)   34,  34,  40 
H  DO  35   r  ■  2»6 
3  5  W (  I )  =  W (  I  )  +  X  (  I  ) 

AX  =  AX  +  X(  1  ) 

GO  TO  30 
36  CALL  CPRNT  (PNORM,  AX,  W»  ERROR,  RONE) 

TF  ( AX  -  Z( 2 ) )   44,  44,  37 
3  7  Z ( 3 )  =  AX 

DO  30   1=1,6 
3  8  S  (  3  ,  I  )  =  W ( I ) 

DL1  =  Z(3  )  -    1(2) 
3°  CONTINUF 

GO  TO  44 

40  PRINT  41 

41  FORMAT  (40H    INFINITE  LOOP  OR  NEEDS  MORE  ITERATION) 
CALL  CPRNT  (PNORM,  AX,  W»  ERROR,  RONE) 

GO  TO  4  4 
4?  PPINT  43 

43  FORMAT  (1QH    CROUT  WENT  CRAZY) 

44  CONTINUE 
CALL  PXIT 
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COMMON  W,  AX,  Bl»  B2»  B3 »  R>  MG  »  A,  C 
END 

SUBROUTINES 


SUBROUTINE  EOR  SOLVING  THE  LINEARIZED  SYST 
OF  EQUATIONS  BY  THE  CROUT  REDOX  METHOD 

LAB^L 

LIST8 

FORTRAN 

SUBROUTINE  CROUT  (A,  C«  X,  w} 

DIMENSION   A(6»6)»  C(6)f  AA(6»6J»  CC(6)i  X( 6  J 

no  ] on   i=i,6 

100  AA(I ill  =  a (  I  ,1  ) 
DO  101   J  =  ?,6 

101  AA( 1 ,J)  =  A( 1  ,J)/A  (  1,1  ) 
DO  102   I  =  2,6 

DO  102   J  =  2,6 

102  AA( I » J)  =  0. 
J  =  2 

103  II  =  J 

DO  1^5   I  =  T  I  ,6 

LI  Ml  =  J  -  1 

DO  104  K    -    1»LIM1 

104  AA(I,J)  =  AA(ItJ)  +  (AA(ItK)  *  AA(KtJ)) 

105  AA(T,J)  =  A(I,J)  -  A  A ( I t  J ) 
IF  ( AA( J, J) )    106,  116,  106 

106  TF  (6  -  J)   110,  110,  107 
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107  I  =  J 

J  =  J  +  1 

JJ  =  J 

DO  109   J  =  JJ,6 

LIM2  =1-1 

DO  108   K  =  l.LIM? 
TOR  AA(ItJ)  =  AA(I»J)  +  (AA(I.K)  *  AA(K.J)) 

AA(  T  *  J)  =  A(  I ♦ J)  -  AA(  I ,J) 

109  AA(I,J)  =  AA(  I  * J)/AA(  I  »  I  ) 
J  =  I  +  1 

GO  TO  103 

110  DO  111   I  =  1*6 
X(  I  )  =  0. 

111  CC( I )  «  0. 

CC(1 )  =  C( 1 ) /AA(  1,1) 
DO  113   I  =  2 » 6 
LIM3  =1-1 
DO  112   K  =  1,LIM3 

112  CC(T)  =  CC(I)  +  (AA(ItK)  *  CC ( K  )  ) 
CC( I )  =  r  < I )  -  CC(  I  ) 

113  CC( I  )  =  CC(  I  )/AA (  T  ,  I  ) 
X(6)  =  CC( 6) 

DO  115   1=1,5 

II  =  6  -  I 

LIM4  =11+1 

DO  114   K  =  LIMA, 6 

114  X (  I  I  )  =  X (  I  I  )  +  (AA(II,K)  *  X(K)) 

115  X(  I  I  )  =  CC(  I  I  )  -  X ( I  I ) 
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SINGULARITY  -  NO  SOLUTION) 


M  =  ] 

no  TO  118 

116  PRINT  1 17 

117  FORMAT  (28H 
M  =  2 

118  RETURN 
END 

SUBROUTINE  FOR  STORAGE  OF  THREE  SETS  OF  VALUES 

OF  AXIAL  LOAD  AND  THE  FOUR  DEFLECTION  COEFFICIENTS. 

*  LA^L 

*  LIST8 

*  FORTRAN 

SUBROUTINE  STORE  (AX,  W,  Z»  S»  LL  ) 

DIMENSION  W(6),  Z(3)»  S(3,6) 

Z(LL)  =  AX 

DO  203   I  =  1.6 
203  S(LL.I )  =  W( I ) 

RETURN 

END 
C  THIS  SUBROUTINE  IS  JUST  A  SHIFTING  PROCEDURE 

C  OF  VALUES  FROM  ONE  STORAGE  LOCATION  TO  ANOTHER. 

*  LABEL 

*  LISTS 

*  FORTRAN 

SUBROUTINE  EXCH  (Z.  5) 
DIMENSION  Z ( 3)  .  S(3,6) 
DO  204   I  =  1,6 
DO  204   J  =  1,2 
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JA  =  J  +  1 
Z04  S( J»  I  )  =  5( JA»  I  ) 
DO  205   IJ  =  1»2 
JB  =  IJ  +  1 

205  Zt IJ)  =  Z( JB) 

DO  2  06   J I  =  1  »6 

206  S(3.JI  )  =  n. 
1(3)     =    0. 
RETURN 

END 

SUBROUTINE  FOR  PRINTING 
LABEL 
LIST8 
FORTRAN 

SUBROUTINE  CPRNT  (PNORM,  AX,  W,  ERROR,  RONE) 
DIMENSION  W(6) 

SUM=W (1 )»*2+W(2)**2+W( 5)**2+9«*(W(3)#*2+W(4)**2 )+25»*W(6)**2 
STRAIN  =  AX  +  3*1416**2  •  SUM/(8.  *  RONE) 
PRINT  200,  PNORM,  AX,  (W<I),  1=1,6),  STRAIN,  ERROR 
200  FORMAT  (F8.2,  F11.5,  6F12.6,  F11.5,  F12.7) 
RETURN 
END 

SUBROUTINE  FOR  COMPUTING  ALL  THE  NECESSARY 

CONSTANTS  IN  THE  MAIN  PROGRAM 
LABEL 
LISTS 
FORTRAN 
SUBROUTINE  CONST  (U,  ALPHA,  °NORM,  R,  Bl,  B2,  B3) 
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DIMENSION  P(13),  R(54)»  Bl(4),  B2(5),  B3(3) 

P(  1  )  =  ALPHA**2 

P  (  2  )  =  1  .  /  P  (  1  ) 

DO  300   I  =  3*7 

Yi  =  (I  -  2)** 2 
3  00  P( I )  =  P( 1 ) / (P ( 1)  +  Yl)**2 

DO  301   J  =  8,11 

Y2  =  (J  -  6)#*2 
301  P(J)  =  P(l)/(1.  +  Y2  *  3(1))** 2 

Pill)    =  P(l)/(4.  *  P(l)  +  9****2 

P(13)  ■  P(l)/(4.  +  9.  *  P(l))**2 

R(l  )   =  P(  1 ) 

R(  2)   =  P( 2 ) 

R(3  )   =  P( 3  ) 

R(4)   =  P(4) 

R(  5 )   ■  P  (  8) 

R(6)   = 

R(7) 

R(8  ) 

R(9) 


P(l )  +       P(2) 
P(l )  +   81.*P(2) 


P(l )  +  625.*P(2) 

■  81.*P(1  )  +  P(2) 

R(10)  =  6?5,*P(1 )  +  P(2) 

R(ll)=  P(l)+  P(8) 

R( 12  )  =  P (1  )  +  81. *P  (  8) 

R( 13 )  =  P(2)  +  P(4) 

P(14)  =  P (2)  +  81.*P(4) 

P ( lc )  =  P ( 1 )  +  3.*P (4) 

R(lft)  =  P(2)  +  3.*P(P) 

R(17)  =  3.*P(1)  +  64.*P(3) 
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R(18  ) 
R(  19) 
R(20) 
R(21  ) 
R(22  ) 
R(23  ) 
R(24) 
R(25  ) 
R(26) 
R(27) 
R(28  ) 
R(29  ) 
R(30) 
R(31  ) 
R(32  ) 
R(33  ) 
R(34) 
R(35) 
R(36) 
R(37  ) 
R(38  ) 
R(39) 
R(40) 
R(41  ) 
R(42  ) 
R(43  ) 
R(44) 
R(45  ) 


3.*P(2)  +  64.*P(3) 

64.*P(°)  +  274.*P(4) 

64.*P(3)  +  274.*P*8) 

P (4)  +  PIP) 

P (3)  +  P(5) 

P  (  3 )  +  P (9  ) 

9.*P ( 1 )  +  64. *P 

9,*P(2)  +  64. *P 

25,*P(1)  +  441. *P 

25i*P(2)  +  44i,*P 

272,*P(3)  +  625. *P 

9,*P(1)  +  64. *P 

9,#P(2)  +  64. *P 

2  5,*P(1)  +  64.  *P 

25.*P(2)  +  64. *P 

64,*P(3)  +  441. #P 

64,*P(3)  +  441, *P 

25,*P(1)  +  576. *P 

25,*P(2)  +  576. *P 

1552, *P(3)  +  28561. *P 

9,#P ( 1  )  +  4.*P  (  2 

4,*PC.)  +  9»*P(2 

36.*P(1)  +  36.*P(2 

?5.*P(15  +  4.*P(2 

4.*P(1)  +  25.*P(2 


36. *P ( 1  )  + 

P(l  )  + 

2.*P( 1 )  + 


P(  2 
36.#P(2 
16«*P ( 3 


3)  4-     25.*P( 

3)  +  25. *P( 
8)  +    256. *P( 

4)  h  256. *P( 
4)  +  625,*P( 
3)  +  1225. *P< 
3)  +  122 5. *P( 
3 )  +  9.*P  ( 
3)  +  9,*P( 
4  )  +  1 6  9  •  *  P  ( 
8)  +  169. *P( 
3  )  +  8  2  8  1  .  *  P  ( 
3  )  +  8  2  3  1  .  *  P  ( 
12)  +  28561. *P( 

+  16.*P(3 )  +  P( 8  ) 
+    16.*P<3)  +         P(4) 

+  225. *P(^ )  +  225. *P( 8 ) 
+  81 ,*P(8 )  +  16.*P( °  ) 
+    81  .  *P(4-  )  +    16.*P(  5  ) 

+  144.*P(9 )  +     9.*P( 13 ) 

+  144«*P(5)  +     9«*P(12) 
+    i+q,*p(^)     +     o.*P(  8  ) 


8) 
4) 
9) 
5) 
8) 
4) 
) 

10) 
6) 
13) 
12) 
13) 
12  ) 
13) 
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R(46  )  =  2.*P(2 )  +  16«*P (3 )  + 
R(47)  ■  16.*P(3)  +  81.*P(4)  + 
R (4fl  )  =   16.*P (3 )  +         P(4)+ 


9.#P(4 )  +    49.*P( 8 ) 

P(8  )  +    81.*P(  10) 
I]  ,*P(  6  )  +    81  »*P(  8  ) 


R(49)  =    4.*P(1)  +   225.*P(2)  +   256.*P(3)  +         P(6) 

R(50)  =  225.*P(1)  +     4»*P(2)  +   256.*P(3)  +        P(10) 

R(51)  =  256.*P(3)  +  2401. *P(4)  +  4096. *P(5)  +  2401.*P(12) 

R(52)  =  256.*P(3)  +  2401. *P(8)  +  4096. *P(9)  +  2401.*P(13) 

R(53)  =   2«*P(H  +  64.*P(3)  +  256.*P(5)  +  225.*P(8)  +  49.*P(12) 

R(54)  =   2.*P(2)  +  64.*P(3)  +  225. *P(4)  +  256.*P(9)  +  49-*P(13) 

DO  302   I  =  1,3 

Y3  =  2  *  I  -  1 

BUI)  =  256.  *  P(l)  *  PNORM/(Y3  *  3»1416**6) 

302  B3(I)  =  -(16.  *  Y3**2 )/( 3. 1416**2 ) 
Rl(4)  =  Bl(2)/3. 

DO  303   J  =  1  ,5 
II  =  2  *  J  +  1 

303  B2(J)  =  4./(3.  *  (1.  -  U**2)  *  P (  I  1  )  ) 
RETURN 

END 

SUBROUTINE  FOR  COMPUTING  THE  COEFFICIENTS 
FOR  THE  USE  OF  THE  CROUT  REDOX  METHOD 

LABEL 

LISTS 

FORTRAN 

SUBROUTINE  SESAME 

COMMON   W,  AX,  Bl,  B2,  R3,  R,  MG ,  A,  C 

DIMENSION       W(6),    "1(4),    B2  (  5  )  ,    83(3),     R(54),     A(6,6),     C(6) 

GO    TO    (4^1  ,     400',     WIG 
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^on    A(l,l)=B3(l)*W(l) 

A(2»l )=B3( 1 )*W(2  ) 

A(3»l )=B3( 2)*W(3) 

A (4, 1 )=B3(2 )*W(4) 

A(5»l )=B3( 1 )*W(5 ) 

A (6» 1 )=B3( 3 )*W(6) 

GO  TO  402 
40  T  A(1  *1 )=B2< 1  )+B3d)*AX+?.*R(6)*W(]  ) ** ?"6 . *R ( 2 ) *W ( 1 ) *W ( 2 ) 

1  -6.*R  (  1  )  *W(  1  )  *W(  3  )+R  (  3  8  )*W  (  2  )  **2+R'.  39  )  *W  (3  )  **2 

2  +^?.MR( 3 )*W( 2 )*W( 3) +  81 .  *R ( 21 )  *W (4)**2-18.*R ( 11 )*W( 2 ) *W ( 4) 

3  -18.*R ( 13 )*W( 3 )*W(4)-6.*R ( 16)*W( 2 )*W( 5  )-6.*R( 15 ) *W  < 3 )*W (6) 

4  +1 6  2.  *R ( 5 )*W ( 4 ) *W( 5 ) +162. *R ( 4) *W (4) *W( 6 ) +R ( 41 ) *W( 5 ) **2 

5  +  R (42 )*W(6 ) **2 

A ( 2, 1 ) =-3.*R ( 2 )*W( 1 ) *#2+2.*R( 3  8 )*W( 1 ) *W  (  2  ) +32.*R( 3 )*W< 1 ) *W ( 3 ) 

1  -18«*R*21  l*WtlJ*W(4)-R(25)*W(g  )**2~2.*R(  24)*W<  2  )  *W  (  -a  ) 

2  +R(40)*W(?)*W(4)-6.*R(16)*W(1)*W(5)+2.*R(18)*W(2)*W(5) 

3  +R(20)*W(3)*W(5)+R(19)*W(3)*W(6)-144.*R(23)*W(4)*W(5) 

4  -9.*R(45)*W(4)*W(6)-R(27)*W(6)**2 

A(3»l )«-3«*R(l )*W(1 )**2+32«*R(3)*W(l)*W(2)+2«*R(39)*W( 1)*W(3) 

1  -18«*R ( 13 )*W< 1 ) *W (4 )-R ( 24) *W(2 ) **2-2«*R ( 25 ) *W ( 2 ) *W { 3 ) 

2  +R  (40)*W(  ?  )*V.' (4)-6#*R(  15  )*W(  1  )  *W(6)+R(  20  )*W(  2  )*W(  5  ) 

?        +R ( 19 )*WI 2  >*W ( 6  >+2«*R ( 17 ) *W I  3 ) *W ( 6 ) ~9»*R (46 )*W ( 4 )*W ( 5 ) 
4        -144«*R(22)*W(4)*W(6)-R(26)*A'(Ei)**2 
A ( 4  » 1 ) =-1 8 . MR (  1 1  ) #W ( 1 J *W< 2 ) -18«*R ( 13 ) *W J 1 J  *W( 3 1 

1  +162  «»R  <  21  )*W  ( 1  )*W  (  4 )  +R  ( 40  )*W(  2  )*  >  I  3 ) +162  »*R<  5 ) *  ■■■■  ( 1  )  *W  (  5  ) 

2  +1  62  .  *R  (  4  )  *W  (  1  )  *W  (  6  )  -1  44  .  *R  (  2  3  )  ■  *l  U )  *  a'  (  5  ) 

-5        -9«*RC45)*W(2)*W(6)-9»*R(46)*W(3)*W(  5)-144.*R(22)*W(3)*W(6) 
4         +1 8 • *R (  1 4 ) *W ( 4 ) *W <  5 ) +1 8 • *R ( 1 2 ) *W ( 4 ) *W ( 6 ) 
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5        +729,*R<;:i)*W(5)*W(6) 
A(5»l)=-6.*R(16)*-W(l)*W(2)+16?.*R(5)*W(i ) #ta (4)+R(l8)*W(2)**2 

1  ■fR(2©)*W<2)*W(3)+9«*R(14)*W(i|.)«#2-144»*Rl23**W4  2  J*W(4) 

2  -o.*R(46)*W(3)*W(4)+2.*R(41)*W(l)*W(5)-2.#R(26)*W(3)*W(5) 

3  +729.*R(21 )*W(4)*W( 6 ) 

A(6»1)=-6.*R(15)*W(1  )*'•/(  3)  +  l62«*R  (  '+)*W(  l  >*W(  4  )+R(  17 )  *WC  3  >*»2 

1  +R (  19)*W(2 )*W( 3 ) +9. *R (  12 )*W(4) **2~V«*R (45) *W( 2  >*W (4) 

2  -144»*R  (  2?  )*W(  3  )*W(  4  )+2.*R  (42  )  *W(  1  )  *W(  6)  -2.*R  <  2  7  )  *W(  2  )  *W  (  6  ) 

3  +729. *R (21 )*W(4)*W( 5 ) 

402  A<l»?)=-3.*R(2)*W(l)**2+2,*R(3fl)*W(l)*W(2J+32«*R(3)*W(l)*W(3) 

1  -18.*R(11)*W(1)*W(4)-2.*R(24)*W(2)*W(3)-R(25)*W(3)**2 

2  +R(40)*W(3)*W(4)-6.*R(16)*W<1)*W(5>+2.*R(18)*W(2)*W(5) 

3  +R ( 20 )*W( 3 ) *W( 5 )+R( 19 )*W( 3 )*W( 6 )-144.*R( 23 ) *W( 4) *W ( 5 ) 

4  -9.*R(45)*W(4)*W(6)-R(27)*W(6)**2 

A( 2  »2  )  =B2 (4  )+B3 ( 1  )*AX  +  R( 3  8 ) *W ( 1  ) **2~2 . *R ( 24 ) * W (  1 )*W( 3 ) 

1  +3«*R(9)*W(2)**2+R(28)*W(3)**?-486«*R(1)*W(2)*W(4) 

2  +R(43)*W«4>*#2+2«*R< 18)*W( 1 )*W( 5J-512»*R(3)*W(3)*W(5) 

3  -2«*R(29)*W(3)*W(6)+R(50)*W(5)**2+R(51)*W(6)**2 

A  (  3  »  2 ) =  1  6  •  *R ( 3 ) *W ( 1 ) **  2-2 . *R ( 24 ) * to ( 1 ) *W ( 2  1-2 • *R ( 25 )*W ( 1 >* w ( 3 ) 

1  +R  (  40  )#W  1  )  *W  (  4  )  +2.  *R  (  73  )  *  W  (  2  )  *W  (  3  )  +R  I  2©  )*W  (  ]  )  *W  (  r-  ) 

2  +R ( 19 )*W( l )*W( 6 )-5l2.*R (3 ) *W ( ?  ) *W( 5 ) -2.*° ( 29) *W(  2 ) *W(6 ) 

3  -2.*R  (  30  )*W(  3  )*W.(  5)-512.-R  (  3  )*W(  3  >*W  (6  )+9.*R(  5  '■  "•  ) 

4  +9.*R(53)*W(4)*W(6)+656.*R(3)#W(5)*W(6) 

A  (  4»  2  )  =-9.*R  (  1  1  )*W  (  1  )  **2+R  (  40  )  *W  (  1  )  *W  (  3  )  -,.43.  *R  (  1  )  *W  (  2  )  #*2 

1  +18»*R(43)*W(2>*W(4)-144.*R(23)*W( 1)*W(5 )-9»*R(31 )*W(5)**2 

2  -Q. *R ( 45 )*W ( 1 ) *W ( 6 ) +9. *R ( 54 ) *W ( 3  )  *W  (  5  )  +  Q  .  *R  (  5  3  )  *W  (  3  )  *W  (  6  ) 

3  -9.#R ( 33 )*W( 5 )*W( 6) 
A(5»2)=-3«*R<16)*W(1 }**2+2«*R< 18 ) * V(  l  ) *W( 2 ) +R ( 20l*W< 1 ) 
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1         -144.*R(23)*W(l)*W(4)-512.*RC*)*W(2)*W(3)-9(33  )*W  (  3  )  **2 
?        +9«*R<54)*W<3)*Wf4)+2.*R(50)*W«2)*W*5)+656»*R<3)*W(3)*W{6) 
3         -18.*R(  31  )*W(  4)*W(  5  )-9.*R(  3  3  )  #  W  ( /+ ) -:.  W(6  )-R(  36)  *W<  6  )**2 

A (6  »2 ) =R  <  1  9 ) *W (  ]  )*W( 3 ) -Q .*R (45 ) *W ( i  )*W ( 4 )-?.*R( ?q ) *W( 2 ) *W ( 3  ) 
1         -256.*R(3)*W(3)**2+9.*R(53)*W(3)*Wf4)-2.*R(27)*A(l)*W(6) 

+2.*R(51)*W(2)*W(6)+656«*R(3)*W(3)*W(5)-0,*R(33) *W (4)*W( 5) 
3        -2«*R (36 ! *W( 5 )*W( 6) 

A(l»3)=-3»*8tl }*WU )**2+2**R(S9)*W<I}*W«3>+32«*R<3)*W<l)*W(2) 

1  -18#*RC13T*W(U*W(4)-R(24)*Wf2  1**2*2 «*R<  25  >*WC  2  )  *W<  3  ) 

2  +R (40 )*W ( 2 ) *W (4)-6,*R (  15 )*W ( 1 )  *W(6  * +  2.*R ( 17 ) *W< 3 )  *W< 6) 

3  +  R  (  20  )*W(  2  )*W  (  5  )  +R(  1  9  )  *W(  2  );;W(  6  )-9.*R  (  46  )*'a;  (  4  )  *W  (  5  )  . 

4  -i44.*R (22 )*W(4!*W(6 )-R (26 ) *W< 5  J**2 
A(2»3)=16«*R(3>*W(1 y**2-2.*R( 25)*Wfl J*W(3»-2«*R(24)*W< 1)*W(2) 

1        +R  (  4  0  )  -*  W  (  1  )  *  W  (  4  )  +  2  .  # R  (  2  8  )  *  W  I  2  J  ■  W  (  3  )  +  R  (  20  }  *W  ( 1 )  *W  (  5  ) 

7  4.R(29)*W!I)*W(6)-2«*R(30)*W(3)*W(5)-512«*R(3)*W(3)*W(6) 

3  -512«*R (3 )*W( 2 ) *W( 5 ) ~2.*R( 29) *W( 2 >*W ( 6  >+9.*R< 54)*W t  4 )  *W ( 5 ) 

4  +9. #R ( 53  )*W ( 4 ) *W ( 6 ) +656. *R ( 3 ) *W ( 5 ) *W ( 6 ) 

A(3  *3  )=B2 ( 2 )+B3 (2 }*AX  +  R ( 39 ) *W ( 1 )**2~2«*R ( 25 ) *W(  1 ) *W« 2 ) 

1  +R ( 28)*W(2 ) **2+3«*R ( 7)*W( 3 )**2~486«*R ( 2 >  *W  ( 3 )*W( 4 ) 

2  +  9.*R (44  )*W ( 4 ) **2+2 .*R ( 1 7 ) *W< \ ) *W{'6 ) -2 t*R< 30) *W  ( 2 ) *W  (  5 ) 

3  -512«*R(3)*W(2 >*W<6 )+R ( 52 )*W<5 )**2*R(49) *W ( 6)**2 

A ( 4  »  3 ) s-9.#R (13) *W ( 1  ) **?+R ( 40 ) *W ( 1 ) *W ( 2 ) -243 . *R ( 2 ) *W ( 3 ) **2 

1  +18«MRf44)*W( 3)*W(4)-9«*R(46)*W( 1 )*W(5 )-9.*R( 32)*W(6)**2 

2  -144«*R(22)*W(  U*W(6y  +  9«*R(54)*W(  2)  *W(5)+9»*R(!        I   (6) 

3  -9»*R(34)*W(5)*W(6) 

A(5*3)=R(2  0)*W(i  )*W<  2)-9**R(46)  *W  (  ]  )#Wf  4  )-236«*      tH  2  )**2 
i        -2«*R  (  30  )*W  (  2  )  *W  (  3  )  +o  .  *R  (54)  *W  (  ?  )  *W  (  4  )  -2«*R  t  26  )  *  '•  (  1  )  *W  (  5  ) 
2        +2»*R(52)*  W  (  3  )  *  W  (  5  )  +  6 r  .         '      I    (6)-9«*R(34  6  ) 


3        -2.*R(35  )*V/(5)*W(6) 
A(6»3)*-3.*R*15)*Wf 1 J**2+2#*R( I7)*wei )*WC3J+R(19)*W< 1)*W(2) 

1  »144«*Rf22)*W( 1)*W(4)-R(29)*W(2)**2-512«*R(3)*W<2)*W(3) 

2  +  9.*R  (  53  )*W  (  ?  )  *W(  4)  +?.*R  (  49  )  *W(  3  )  *W(  6  )  +656.*R  (  3  )  •  '■  (  '-M  ■  '  I  5) 

3  -9  .  *R  (  3  4  )  «W  (  4  )  *W  (  5  )  -  1  P  .  *R  (  3  2  1  *  W  (  4  )  *W  (  6  )  -R  (  3  5  )  *  -■  (  5  )  *  l 
A ( 1  i 4  )  =  1 6  2 . MR ( 2 1 ) *W (  1  ) *W ( 4 ) - 1 R . *R ( 1 1 ) * W ( 1  ) * W ( 2 ) 

1  -1 8.*R( 13 )*W( 1 ) *W<3 )+R ( 40) *W ( 2 t*W< 3 ) +162 «*R < 5 )*W (1)*W<5) 

2  +  162.*R  (4)*W(  1  )  *W(6  )  -144.*R  (  2"M  *a<(  2  )  *W  (  5  ) 

^        -9.*R(45)*W(2)*W(6)-9«*R(46)*W(3 )*W< 5 ) -144 »*R ( 22 ) *W ( 3 ) *W ( 6 ) 

4  +  18«*R(14)*W(4)*W(5)+18«*R(12)**/(4)*W(6) 

5  -729. *R (21  )*W( 5 )*W( 6 ) 

A ( ?,4)=-9.*R ( 1 1 )*W( 1 ) **2+R (40 ) *W< 1 ) *W( 3 )-243«*R( 1 )*W( 2 ) **2 

1  +2*tR(43)*W«2>*W(4)-144»*R(23)*W(l)*W( 5 ) -9«*R (45 ) *W( 1 )*W< 6 ) 

2  +9.*R(54)*W(3)*W(5)+9#*R(53)*W(3 ) *W ( 6 ) -9«*R < 3 1 ) *W ( 5 }**2 

3  -9.*R ( 3  3 **WI 5 )*W( 6) 

A(  3»4  )  =-9.*R  (  1  3  )*W  (  l  )  **2+R  (  40  )  *W(  l  )*W(  2  )-?43«*R  (2  )*Wt3)*  *2 

1  +18.*R ( 44)*W ( 3 ) *W( 4 ) -9.*R ( 46)*W ( l )*W ( 5 ) -9 «*R ( 34 ) *W ( 5 ) *W ( 6 ) 

2  -144.*R  (22  )*W(  1  J*W(f>  )+o,*R  (  54)  *W(  2)  *W  (  5  )  +9.*R(  53  )  *W  <  2  )  *  a'  (  6  ) 

3  -9.*R (32 )*W(6 )**2 

A(4»4)=81  .  *B2  (  1  )  +  R3  (  ,?  )*AX+81.#R  (  2  1  )  *W  (  1  )  **2  +  9  «*R  (  43  )  *W  (  2  )  **2 

1  +243«*R<61*Wt4)**2+9«*R(44)*W( 3 )**2+18»*R( 14) *W ( 1 ) *W( 5 ) 

2  +1 8 • MR  « 1 2 ) *W  C 1 >*W 4 6 ) +8 1 t*R ( 47 ) *W ( 5 ) **2  +  8 1 • *R ( 48 ) *W ( 6 ) **2 

3  +2592«*R(3)*W(5)*W(6) 

A(  5  »4)=R1  .*R  (  5  )*W(  1  )#*2+l8**R  (  14)  *W  I  1  )  *  W  < /•  )  +Q.*R(  54)  *W  (  ?  )  *W  (  3  ) 

1  -144.#R (23 )*W(  1 )*W( 2 )-9.*R (46) *W( 1 ) *M  3 ) 

2  +729«*R(21 )*W( 1 )*W( 6 )-18.*R ( 31 ) *W(2 ) *W ( 5 ) 

3  -9.*Rt33)*W(2)*W(6)-9«*R(34)*W(3 )*W( 6 ) +162 **R ( 47 ) *W( 4)*W ( 5 ) 

4  +2592«*R(3)*W(4)*W(6) 
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A(6*4)=81.*R(4)*WM  )  **?  +  lP.*R  M?  )*W(  l  J*W(4)-9«*R(45)*W«1 )*W<2) 

1  -144.*R(22)*W(1)*W(3)+9.*R(53) *W ( 2) *W( 3 ) 

2  +729.*R(2I)*Wt  1  )*W(5 J*9«*Rt33)  *W  (  2  J  *  W  (  5  )  -9  .  *R  (  3  4  )  *W  (  3  )  *W  (  5  ) 

3  -18»*M32)*W(3)*W(6 )+2  592  **R(3) *W<4 ) »W ( 5 ) 

4  +162.*P(48)*W(4)*W(6) 
AM»5)=-6.*Rn6)*W(l)*W(2)+16?.*R(5)*W(l)*W(4)  +  ?.*R(41)*W(l  )  *W  (  5) 

1  +R  (  18  )*W(  2  )  **2  +  R(  2  0  )  *W(  2  )  *W  (  3)-144«*R  (  2  3  )*W<  2  )*W(  4) 

2  -9.*R(46)*W(3)*W(4)+9.*R(14)*W(4)**2-2.*R(26)*W(3)*W(5) 

3  +729.*R(21 )*W(4)*W(6 ) 

A  (  2  >5  )=-3.*R ( 16  )*W( 1  )**2+2 .*R  *  18 ) *W( i ) *W< 2 ) +R ( 20)*W( l )#W ( 3 ) 
1        -144.*R(23)*W(1)*W(4)-RC30)*W(3)**2~512«*R<3)«W(2 )*W(3) 
?         +  9.#R ( 54 )*W ( 3 ) *W( 4) +  ?.*R ( 50 ) *W( 2 ) *W ( 5 ) +6  5  6.*R ( 3 ) *W ( 3) *W  ( 6 ) 
3         -18«*R( 31 )*W( 4) *W( 5 )-9,*R ( 33)*W (4)*W ( 6 )-R( 36)*W(6 ) **2 

A  (  3  »  5  )  =R  (  2  0  )  *W  (  1  )  *W  (  ?  )  -Q  •  #R  (  46  )  *W  (  1  )  *W  (  4  )  -2  •  *R  (  26  )  *W  (  1  )  *W  (  5  ) 

1  -2  56#*R  (3  )*W(  2  )  **2-2  #*R  (30  )*W(  2  J  *W  (  3  )+9.*R(  54)  *W  (  2  )*W(4  ) 

2  +  2.#R ( 52  )*W( 3 )*W( 5)+656.*R ( 3 )  #W  ( 2 )*W< 6 >-9»*R( 34) *W(4)*W ( 6) 

3  -2.*R(35 )*W(5)*W( 6) 

A ( 4  » 5 ) =  8 l . * R ( 5  )  *W ( i  ) **2-l 44«*R ( 2 3 ) *W ( ]  ) *W ( 2 ) ~9 • *R ( 46 ) *W ( ]  ) *W ( 3 ) 

1  +18«*R{I4)*W< l)*W(4)+?29«*R(2iy*W(iy#W(6) 

2  +9.*R(54)*W(2)*W( 3)-18.*R « 31)*WC2 )*W< 5  )-9.*R< 33 ) *W ( 2 ) *W ( 6 ) 

3  -9.*R (34)^W( 3 )*W(6)+162«*R (47) *W(4) *W  ( 5) 

4  +2592.*R(3)*W(4)*W(6) 

A( 5*5 ) =B2 ( 5 )+B3 ( 1 )*AX+R( 41 )  *W  ( 1 >**2~2«*R ( 26 )*W( 1 ) *to( 3 ) 

1  +R  (  5  0  )*W(  2  )  #*2+D(  5  2  )*W<  3  )  **2~18#*R(  31  )  *W  (  2  )  *W  (  4  ) 

2  +8  1  ,*R(  47)  *W{  4)-** 2 -2  .*R  (35  J*W(  3 ) *W t 6  )+3.*R ( 10)  *W  (  5  )  **2 

3  +  P  (  37  )*W(6 )**2 

A (6»5 ) =72°#*R( 21 )  *W ( 1 )  *W  <  4 ) +6  56 «*R < 3 ) *W < 2 ) *W( 3) 
1        -9«*R(33)*W(2)*W(4)-9«*R(34)*W(3)*W(4)+1296«*R(3)   (4)* 
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2        -2.»R(35)*Wf3)*W«5)-2«*R<36)*W<2)*W<6)+2»*R(37H       ( 5 ) 
A( l*6)*-6«*R(15)*W(l}*W(3)+162**R(4>*W(l )*W<4)+2»*R(42)*W(l )*W(6) 

1  +R  ( 17  )*W  {  3  J**2+R  ( 19  J  *W  t  2  J  *W  13  )  -§»*R  (  45  )*W  (  2  ) *W  (  4  ) 

2  -144*»R(22>*W(3)*W<4)+9#*R(12)*W(4***2-2.*R(27)*W(2)*W<6) 

3  +729.*R(?1)*W(4)*W(5) 

A  (  2  »6  )  =R  (  19  )*W  (  1  )*W<  3  )-Q.*R  (45  )  *W  (  1  )*W(  4)-?#*R(  ?7)  *W<  1  >*W  (  6  ) 
1        -256«*R<3)*W(3)**2-2«*R(29)*W( 2)*Wi  3 )+9«*R( 53)*W(3)*W(4 ) 
?        +2»*R<51 >*W(2)*WJ6)*656»*Rf 3)*W(3)*W(5 )-9«*R( 33)*W (4)*to 
3        -2.*R (36 )*W( 5 )*W( 6) 

A( 3  »6)  =-3.*R ( 15  )*W( 1 )*#2+R ( 19) *W( 1 )*W( 2 )+2.*R ( 17) *tf  ( 1 )*W{ 3 ) 

1  -144#*R(2?)*W(1)*'W<4)-R(29)*W(2)**2-512.*R(3)*W(2)*W(3) 

2  +9.»R!53  >*W(2)*W(4)+2«*RJ4t>*W( 3)*W(6)+656«*R(3)*W(2)*W( 5) 

3  -1 8 . *R ( 3? ) * W ( 4 ) *W ( 6 ) -9 . *R ( 3  \) #W I 4 ) *W  <  5  J -R ( 35 ) *W ( 5 ) **2 
A(4»6)«81«*R<4r*W(l )**2»9«*R(45)*W(1)*W(2)-144»*R(22)*W( l)*W(3 ) 

1  +  1 8 •  *R  (  12  ) *W  {  1  )'*W  ( 4 )  +729  **R  (  21 )  *W  ( 1 J  *W  (  5  J 

2  +9.*R (53  )*W(2)*W( 3)-18.*R( 32)*W(3 )*W(6)-9.*R( 33)*W(2)*W( 5) 

3  -9.*R(34)*W(3)*W(5)+162«*R(4fl)*W(4)*W(6) 

4  +2  592.*R (3 )*W( 4)*W( 5 ) 
A(5»6)=729.*R(21)*W(1)*W(4)+656.*R(3)*W(2)*W (3) 

1  -9#*R(33)*W(2)*W(4)-9.*R( 34)*W(3 )*W( 4) +1296«*R ( 3 )*W(4)**2 

2  ■*2.*R(36)*W(2)*W<6)-2«*R(35)#W(3  )*W<  5  )  +  2.*R  (  37  )*W  (  5  )  *W(  6  ) 
A  (6»6  )  =B2  I  3  )+B3  (3  )  *AX  +  R(  42  l*W  (  1  )**?.-?_  .*R  (27  )*W<  1  )  *W<  2  ) 

1  +R  (  5  1  )  *  W  (  2  )  **  2  +  R  (  4Q  )  *W  (  3  )  *  *  ?  - 1  8  •  *  R  <  3  ?  )  *  W  (  3  )  *  W  (  4  ) 

2  +81«*R<48)*W(4)**2-2«*R<36)*Wt  2)*W«5)+R( 37)*W( 5 

3  +3.*R{8>*W(6)**2 

01 =P1 ( i ) -R2 ( 1 )*W( 1 ) -83 ( 1 )#AX*W (1)-R(6 >*W( 1 ) **3 

1  +3«*R(2 )*W(1 J**2*W(2)+3t*R<l)*W(l )**2*W(3)-R(38)*W( 1 )*W(2)**2 

2  -R(39)*W(l)*W(?)**2-32.*R(3)*W(l)*vV(2)*W(3) 
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3  -81.*R ( 2  1 )*W( 1 )*W(4)**2+18«*R ( 11 )*W< 1 )*W ( 2 )*W< u ) 

4  +13«*Rll3»»W(l?*Wf  3  )*Wf 4)+R<24)*Wf 2)**2*W(3) 

5  +R(25)*Wt2J*W(3) **2-R <40 >*W (2 ) *W ( 3 ) *W<4) 

6  +  6  .  *R  (  1  6  )  *W  ( 1  )  «  'a  (  ?  )  * w  (  5  )  +f  •  *R  (  1  3 )  #  .v  ( 1 )  *W  t  3  )«  W  (  6  ) 

7  - 1 6  ? . *R ( 5 ) * W ( 1  ) # W ( 4 ) * W ( 5 ) - 1 62 • *R ( 4  )  *W  (  1  )  #W  {  4  )  *W  (  6  ) 

El  s-R  (  41.  )  *W  (  l  )  *W  (  ^  )  **2-R  (  42  )  *W  (  l  )  *W  (  f  )  *  *?-R  (  }  fl  )  *W  (  ?  )  *#2*W  (  «j  ) 

1  -R 1 17)*W I S***2*W( 6 ) -R ( 20  >*W (2 ) *W ( 3 ) *W ( 5 ) 

2  -R(19)*WC2>*W(3>»Wt6)+144«*R(23>*W(2 )*WC4)*W(5 ) 

3  +o.*R  (  45  )  *W  (  2  )  *W(  4  )  *W  (  6  )  +  9.#R  (  45)  *W  (  3  )*>.'  (  4  )  *W  (  5  ) 

4  +144.*R  (  22  )  *W(  3  )  *W(4)  *W  (  6)-9»*R(  14)*W(4)  **2*W(  5  ) 

5  -9«*R(  12)*W(/- )**2*W(6)+R(27)*W(2)*W(6)**2 

6  +R  (  26  )*W(  3)*W(  5  )  **2-729  .  #R  (  2  1  )  *W  (  4  )  *W  (  5  )  *W  (  6  ) 
D2=R1  (2 ) -3  2 ( 4)*W( 2 ) -B3( 1  )*AX*W(? ) +R  <  ? ) *W< 1 ) **3 

1  -R(38)*W(1)**2*W(2)-16,*R(3)*W(1)**2*W(3) 

2  +9.*R(HD*W(l)**2*W(4)+R(25)*W(l)#W(3)-»-*2 

^  +2 • *R ( 24 ) *W ( 1  ) *W ( 2 ) *W ( 3 ) -R ( 40 ) *W ( 1  ) #W ( 3 ) »W ( 4 ) 

4  -R(9 ) *W< 2 }**3-R( 28 ) *W  (  2 ) *W ( 3 )**2-R( 43 )  *W ( 2 )*W(4 )**2 

5  +24?«*R(1)*WC?)**?*W(4)+3#*R(16 ) *W < 1 ) **2*W ( 5 ) 

6  -2«*R{]  8D*W(  1  )*W(2)*W<5)^-R(20J*  V(l  > *W  ( 3  >  *W  (  5  ) 

7  -R( 19)*W( 1 )*W( 3 ) *W( 6)+144.*R( 23)*W( 1 )*W( 4) *W  (  5 ) 
F?  =  9  ,*R  (  4-5  )  #W  (  1  )*W  <  4  )  *W  (  6  )  +R  (  27  )  *W  (  l  )  *W  (  6  )  **2 

1  +R( 30 )*W ( 3 ) *#2*W( 5 ) +2  56.#R ( 3 }*W ( 3 ) **2*W( 6 ) 

2  +  512»*R(3  )*Wf  2  J*Wf  3)*W(5)  +  2«*R(29)*W<2)*W(  3)r*W(61 

3  -9«*R ( 54 )*W(3 )*W( 4 ) *W ( 5 ) -o . *R ( 53 ) *W ( 3 ) *W ( 4 ) *W  (  6 ) 

4  -Rf 50)*WC2? *W<5 )**2-R <  5  1 >*W(2 ) *W(6 )**2 

5  -656#*R( 3}*W(3 J*W(5 )*W(6)+9«*R(31)*W(4)*W( 5)**2 

6  +9 •  *R  (  33  \ *W  ( 4 )  *W  (  5  )  *W  (  6  )  +R  ( -36  )  *  W  (  c  )  *  '-'  ( f. )  ** 2 
D3«B1  (  2  )  -P.  2  (  2  )  *W(  ?  )  -B3  <  2  >*AX*W(3  )  +R  (  1  )  *W  (  1  )  #*3 
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1  -16.*R ( 3 ) *W  U )**2*W< 2 )-R ( 39 ) *W (1 )* *2*W I  3 ) 

2  +  o  .  *  R  (  i  3  )  *  W  ( 1  )  *  *  2  *  W  (  4  )  +  R  (  2  4  )  #  a1  (  l  ) #  W ( 2 ) **2 

3  +2 .*R ( 2  5 ) *W ( l  ) *W ( 2 ) *W ( 3 ) -R ( 40 ) *X ( 1  ) *W ( 2 ) *W ( 4 ) 

4  _r ( 2  8 )*W ( 2 ) **2*W ( 3 )-R ( 7 ) *W< 3 )**3+243«*R< 2 ) *W ( 3 ) **2*to '  ( 4 ) 

5  -9.*R( 44)*W (3 )*W(4)**2+3«*R ( 15  >*W< 1 ) **2*W( 6 ) 

6  _r ( 2  0 ) * W ( 1  ) * W ( ? ) * W ( 5 ) -R (  1  9 ) * W (  1 ) * W ( 2 ) *W ( 6 ) 

7  -2.*R ( 17 ) *W ( 1 )*W< 3 ) *W(6 ) +9.*R ( 46)*W ( 1 )*W ( 4 )*W( 5 ) 

E3  =  ,44.*R(22)*W(l)*W(4)*W(6)+R(26)*W(l)*W(c!;)**2 

1  +256.*R(3)*W(2)**2*W(5)+R(29)*W(2  )**2*W<  6) 

2  +2,*R(30)*W(?)*W(3)*W(5)+512«*R<3)*W(2)*W(3)*W(6) 
■3  -9,#R  (  5  4  )*W  (  2  )*W  (  4)  *W  (  5  )  -o#*R  (  5  3  )*W(  2  )*W  (  4  )*W<  A  ) 

4  _R( 52 )*W ( 3 )*W( 5 )**2-R (  49  )  *W ( 3 ) *W ( 6 »**2 

5  -656.*R(3)*W(2)*W(5)*'"M6)+9#*R(32)*W(4)*W(6)**2 

6  +q.*R ( 34)*W(4)*W( 5)*W(6 )+R ( 35 ) *W ( 5 ) **2*W ( 6 ) 

D4  =  R1  (4  )-81 .*B2( 1 )*W(4)-B3 I  2 ) *AX*W ( 4 ) +9 . *R (  1  1 )*W( 1 )**2*W(2 ) 

1  +9.*R ( 1 3 ) *W ( 1 ) **f*W ( 3 ) -8 1 .*R ( 2 1 ) *W ( 1 ) **2*W ( 4 ) 

2  -R(40 )*W ( 1 )*¥{ 2 ) *W( 3 )+81 «*R (1)*W(2 ) **3+81 »*R(2)*W(3) **3 

3  -^,*Re43)*W«2j**2*W(4)-8i«*R(6)*W(4)**3-9«*RU4)*W(3)**2*W{  4) 

4  -81.*R  ( 5 ) *W ( 1  )**2*W<  5  )-81«*RI4  >*W(  1  )*#2*W  (  6) 

5  +144.*R ( 23 ) *W ( 1) *W( 2 ) *W ( 5 ) +Q.*R( 4* ) *W ( 1 ) *W ( 2  )*W (6  ) 

6  +9,*Rf 46)*WI 1 )*Wf 3)*W<5 )+144«*Rt 22 )*W( 1) *W(3)*W(6 ) 
F4  =  - 1  8 . #R (  1 4 ) *W ( 1  ) *W ( 4 ) # W (  5 ) -  1  8 . *R (  1 2 ) *W ( 1 ) *W ( 4 ) #W ( 6  ) 

1  -7?q  ,*R  (  2  1  )*W  ( ]  )  *W  (  5  )  *W  t  6)  ^9»*R  <  54)  *W  (  2  )  *W  (  3  )*W  (  5  ) 

2  -o.*R( 53 )*W (2  )*W( 3 ) *W ( 6 )+9.*R ( 31 }*W( 2 )*W ( 5 )**2 

3  +9,*R(32 )*W(3 )*W(6)**2+9«*R(33 ) +W(2)*W( 5 )*W(6) 

4  +9.*R ( 34) *W  f  3  ) *W ( 5 ) *W ( 5 ) -8 1 ,*R ( 47 ) *W ( 4 ) *W ( 5 ) **2 

5  -8 1 • *R ( 48 ) * W ( 4 ) * W  <  6 ) **2- 2  5 9 2 • *R ( 3 ) *W ( 4 ) *W ( 5 ) *W ( 6 ) 
05=61 ( 3 )-B2 ( 5)*W(5 ) -B3 ( 1 )*AX*W(5 )+3.*R( 16)*W( 1)**2*W(2 ) 
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1  -81  .*R  (  5  )  *W(  l  )**2*W(4  )-R  (  1  8  )*W  (  1  )*W(  2  )**2-R(  20)  *W(  1  )  *W  (  2  )  *W  (  3  ) 

2  -9 . *R ( 1 4 ) *W ( 1  ) *W ( 4 ) **2  +  1 44 . *R ( 2  3 ) *  W ( 1  ) *W ( 2 ) *W ( 4 ) 

3  +  9.*R ( 46 )*W( 1  )*W(3)*W(4)+256»*R< 3) *W(2)**2*W (3 ) 

4  +R(30)*W(2)*\.'(3  )  **2-9»*R  (  54  )  *W  (  2  )  *W(  3  )*W  (  4  ) 

5  - R ( 4 1  ) *W ( 1 ) **2*W ( 5 ) +2  »*R ( 2 6 ) *W ( 1  > * W ( 3 ) *W ( 5 ) 

6  -729«*R( 21)*W(1)*W(4)*W(6)-R( 50) *W ( 2 )**2*W ( 5 ) 
f5  =  -R ( 52)*W ( 3 ) *#2 *W ( 5 ) -656 . *R ( 3 ) *W ( 2 ) *W ( 3  I *W ( 6 ) 

1  +18«*R<31)*W(2)*W(4)*>«'(5)+9,*R(33)*W(2)*W(4)*W(6) 

2  +  9,*R(34)*W(3)*W(4)*W(6)-81»*R(  '+7)  *W(4)**2*W(5  ) 

3  -1296.*R( 3)*W(4) **2*W(6 ) +R ( 36)*W< 2  )*W(&) **2 

4  +2**R(35 )*W(3)*W(5)*W(6 )-R( 10)*W( 5 ) **3-R ( 37 ) *W ( 5 ) *W( 6 ) **2 
D6=B1 (3)-B2(3)*W(6)-B3(3)*AX*W(6)+3«*R( 15)*W( 1)**2*W(3) 

j    -81«*R(4)*W(1 )**2*W(4)-R(17)*W(1 )*W(3)**2-R< 19)*W( 1)*W(2)*W<3) 

2  -q,*R (12 )*W(1 )*W( 4)**2+9«*R (45 )*W( 1 )*W( 2 )*W< 4) 

3  +i  44«*R ( 2? ) *W( 1 )*W( 3)*W (4)+R ( 29) *W(2 )**2*W ( 3  ) 

4  +256. *R( 3)*W(2)*W(3)**2-9«*R< 53)*W( 2 )*W( 3  »*W(4J 

5  „r ( 42 ) *W ( 1 ) #*2*W ( 6 ) +2 • *R ( 2  7 ) *W ( 1 ) * W ( 2 ) *W ( 6 ) 

6  -729  «*R ( 2 1 ) *W( 1 ) *W ( 4 ) *W i 5 ) -R ( 5 1 ) *W ( 2 ) **2*W ( 6 ) 
E6  =  _r ( 49 ) #w ( 3 ) **2*W ( 6 ) -656  «*R ( 3 ) *W ( 2 ) *W ( 3 ) *W ( 5 ) 

I  +Q,#R(33)*W(?)*W(4) *W< 5 )+9«*R( 34)*W(3)*W(4)*W( 5 ) 

?         +18«*R(32)*W(3)*W(4)*W(6)-1296«*R(3 )*W( 4) **2*W( 5 ) 

3  -81«*R(48)*W(4)**2*W(6)+R(35)*W(3)*W(5)**2 

4  +?  #*R  (  36  )*W  (?  )  *W  (  5)*W(6  )-R  (  37)*W(  5  )  **2*W  (  6  )  -R  (  -p  )  *W  (  6  )  **3 
C(l )  =  Dl  +  El 

C(2)  =  02  +  E2 
C ( 3 )  =  0  3  +  F  3 
C(4)  =  D4  +  F4 
C  (  5  )  =  Dri    + 
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C(6 )     =    06    +    E6 
RETURN 

END 
DATA 
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